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Abstract 


We review the current status and recent progress of microscopic many-body approaches and 
phenomenological models, which are employed to construct the equation of state of neutron stars. 
'The equation of state is relevant for the description of their structure and dynamical properties, 
and it rules also the dynamics of core-collapse supernovae and binary neutron star mergers. We 
describe neutron star matter assuming that the main degrees of freedom are nucleons and hyper- 
ons, disregarding the appearance of quark matter. We compare the theoretical predictions of the 
different equation-of-state models with the currently available data coming from both terrestrial 
laboratory experiments and recent astrophysical observations. We also analyse the importance of 
the nuclear strong interaction and equation of state for the cooling properties of neutron stars. We 
discuss the main open challenges in the description of the equation of state, mainly focusing on the 
limits of the different many-body techniques, the so-called *hyperon puzzle," and the dependence 
of the direct URCA processes on the equation of state. 
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1 Introduction 


Neutron stars (NS) harbor unique conditions that challenge the physical theories of dense matter [1]. 
Beneath a thin stellar atmosphere, a NS interior consists of three main layers, i.e. an outer crust, an 
inner crust, and a core, each one characterized by different physical conditions. At densities larger 
than 1.5 x 1014 g/ cm’, nuclear matter forms a homogeneous liquid composed of neutrons plus a certain 
fraction of protons, electrons, and muons that maintain the system in beta-equilibrium, whereas in the 
deep interior, at still higher densities, strange baryons and even deconfined quarks may appear. On 
the other hand, moving from the core to the exterior, density and pressure decrease. When the density 
becomes lower than approximately 1.5 x 10% g/ cm’, clusterized nuclear matter sets in. The positive 
charges concentrate in individual clusters of charge Z and form a solid lattice to minimize the Coulomb 
repulsion among them. The lattice is embedded in a gas of free neutrons and a background of electrons 
such that charge neutrality is maintained [1, 2]. This layer is the so-called inner crust, and the nuclear 
structures may adopt non-spherical shapes (usually termed “nuclear pasta”) in order to minimize their 
energy [3-5]. At lower densities, neutrons are finally confined within the nuclear clusters and matter is 
made of a lattice of neutron-rich nuclei permeated by a degenerate electron gas. This layer is known as 
the outer crust [3] and extends from the neutron drip density Parip © 4 x 10% g/cm” to about 104 g/cm” 
at the surface. 

It is usually believed that asymmetric nuclear matter forms the interior bulk part of NSs. Despite 
infinite nuclear matter being obviously an idealized physical system, the theoretical determination of the 
corresponding equation of state (EoS) is an essential step towards the understanding of the NS physical 
properties, and the comparison of the theoretical predictions with the experimental observations can 
provide serious constraints on the nuclear EoS. Unfortunately, only indirect observations of NS structure, 
including their masses and radii, are possible. However, the astrophysics of NSs is a rapidly developing 
research field, in view of the observations coming from the currently operating satellites (NICER, 
Neutron star Interior Composition Explorer, [6]) and the gravitational-wave (GW) laser interferometers 
(AdvLIGO [7, 8], and Virgo [9]), and it is now possible to confront the theoretical predictions with 
more and more stringent phenomenological data. 

Heavy ion reactions is another field of research where the nuclear EoS is a relevant issue. In this case, 
the difficulty of extracting the EoS is due to the complexity of the processes, since the interpretation of 
the data is necessarily linked to the analysis of the reaction mechanism. An enormous amount of work 
has been done in the last two decades in this field, and some indications about the main characteristics 
of the EoS emerged, but limited to a density range slightly above the nuclear saturation density, and 
therefore far away from the high density range typical of the NS interior. However, a direct link 
between the two fields of research can also be impeded by the fact that the typical time scale of heavy 
ion reactions is enormously different from the typical NS time scale. Indeed, nuclear matter inside NSs 
is completely catalized, i.e. it is quite close to the ground state, reachable also by weak processes. In 
heavy ion reactions the rapidity of the dynamical evolution can hinder the weak processes which relax 
the system towards such a catalized state, and therefore the tested EoS can differ from the NS one, 
especially at high density. 

Further astrophysical scenarios characterized by high nuclear density and large temperature, and 
where the EoS plays a relevant role are i) core-collapse supernova explosions (CCSN), which lead to 
the formation of a very hot proto-neutron star (PNS) and subsequently to a cold NS or to a black hole 
(BH) [10-13], ii) merging of NS in close binary systems, NS-NS and NS-BH [14-16]. The dynamical 
evolution and the structure of the forming final compact object are strongly determined by the nuclear- 
matter EoS, as well as the nucleosynthesis conditions and the emitted neutrino spectra. Therefore there 
is a deep connection between the gross features of these astrophysical phenomena and the underlying 
microphysics ingredients which arise from the interactions among the particles. 

Within this review we will consider only nucleonic and hyperonic degrees of freedom, and disregard 
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the appearance of quark matter, since this is a very complex topic which deserves a review by itself; we 
refer the interested reader to Ref. [17, 18]. 


The main goal of this work is to review current theoretical approaches for the description of NS 
matter that can yield EoSs relevant for compact stars. After an extensive reexamination of the baryon- 
baryon interaction in Sec. 2, where all the different theoretical approaches are discussed, we concentrate 
in Sec. 3 on the different theoretical methods currently used for the construction of the EoS. We illustrate 
both ab-initio and phenomenological models, at zero and finite temperature, and confront them with 
experiments. After a short review of the crust EoS in Sec. 4, we will illustrate in Sec. 5 the role of 
hyperonic degrees of freedom in nuclear matter, and its consequences on the NS structure. In particular, 
the so-called “hyperon puzzle” is reviewed. Sec. 6 is devoted to illustrate the current nuclear physics and 
astrophysical constraints available so far, and show how the different NS models behave with respect 
to the laboratory and observational data. Finally, in Sec. 7 we discuss how the EoS determines the 
cooling behaviour of NSs, emphasizing the roles of superfluidity and nuclear pairing gaps. Conclusions 
are drawn in Sec. 8. 


2 Baryon-baryon interaction 


The features of the baryon-baryon interaction [nucleon-nucleon (NN), hyperon-nucleon (YN), hyperon- 
hyperon (YY)], in particular the presence of a hard repulsive core, strongly determine the behaviour 
of the nuclear medium, and are the basic ingredient of calculations performed within microscopic ap- 
proaches. The dominant component is the two-body interaction, but higher-order forces can be im- 
portant, e.g. the three-body force (TBF) which is crucial in reproducing the saturation properties of 
nuclear matter. The first model of the strong force between hadrons was proposed by Yukawa in 1935 
[19], and this model was used to describe the nuclear force between nucleons mediated by the 7-meson, 
whose mass was related to the range of interaction, thus restricting the potential to a finite range, which 
is the essential point. Later on, with the discovery that the fundamental theory of strong interactions is 
quantum chromodynamics (QCD) and not meson theory, all attempts to derive a proper theory of the 
nuclear force had to be formulated again. Soon it became clear that the nuclear Hamiltonian should in 
principle be derived from QCD, but this is still a very difficult task to solve, and therefore we have to 
resort to alternative approaches. Around 1990, a major breakthrough occurred when the concept of an 
Effective Field Theory (EFT) was applied to low-energy QCD [20, 21]. This scheme is also known as 
Chiral Perturbation Theory (ChPT) and allows to calculate the various terms that make up the baryon 
force systematically order by order, thus generating not only the force between two baryons, but also 
many-baryon forces, all on the same footing. Nowadays, the bare baryon interactions can be basically 
divided into four classes: i) phenomenological interactions; ii) chiral effective field theories (xEFT); iii) 
resonating group method (RGM) models that include explicitly the quark-gluon degrees of freedom, 
and iv) renormalization group methods. A major progress has been also achieved recently on derivation 
of the baryon-baryon interaction from lattice QCD calculations. We will briefly discuss all of them in 
the following. 


2.1 Phenomenological models 


Theoretical approaches to construct phenomenological forces include meson-exchange models and po- 
tential models. In both cases, quark degrees of freedom are not treated explicitly but are replaced by 
hadrons — baryons mesons and their resonances — in which quarks are confined. 


2.1.1 Meson-exchange models 


Meson-exchange models are based on the Yukawa idea according to which the strong interaction among 
the different baryons is assumed to be mediated by the exchange of mesons between the baryons. The 
long-range part of the interaction is mediated by the exchange of pseudoscalar mesons (7, K, m, 7), 
whereas scalar mesons (0, k, 0) describe the intermediate-range attractive part, and the vector mesons 
(p, K*, w, $) mediate the short-range repulsive contribution. The various models differ mainly in the 
mesonic content and the treatment of two-meson exchange contributions, but they are very successful in 
describing NN scattering data and phase shifts at laboratory energies S 350 MeV, and also the deuteron 
properties, even if discrepancies between the results of different groups do still exist [22]. Very refined 
and complete phenomenological models have been constructed for the NN interactions, e.g. the Paris 
potential [23], the Bonn potential [22, 24], and the Nijmegen potentials [25, 26]. Meson-exchange theory 
has been employed by the Nijmegen [27-31] and Júlich [32, 33] groups to construct also YN and YY 
interactions. 

Following symmetry principles, simplicity and physical intuition, the most general interaction La- 
grangians that couple the meson fields to the baryon ones can be written as 


Ls = gll. (1) 
Los = gps Vin Vo” , (2) 
Ly = g Vy Vo + goy (9,9 — 0,90) , (3) 


for scalar, pseudoscalar and vector coupling, respectively. Alternatively, for the pseudoscalar field one 
can write also the so-called pseudovector (pv) or gradient coupling, which is suggested as an effective 
coupling by chiral symmetry 

Lyo = gp, Y y? y à, BP) : (4) 


In the above expressions V denotes the baryon fields for spin-1/2 baryons, ll, dipl and PT are 
the corresponding scalar, pseudoscalar and vector fields, and the g’s are the corresponding coupling 
constants that must be constrained, when possible, by scattering data. Note that the above Lagrangians 
are for isoscalar mesons. To obtain the Lagrangians describing the couplings with the isovector mesons, 
the fields ® should be simply replaced by 7 - ® in the previous expressions, with 7 being the usual 
isospin Pauli matrices. 

A typical contribution to the baryon-baryon scattering amplitude arising from the exchange of a 
certain meson ® is given by 


ü(pi)gaViu(pi) Peti(p5) gel 2u(pa) (5) 
(p — 24)? — må | 


(pips|Va|pipz) = 


where På /[(p, — p1)? — m2] represents the meson propagator, my is the mass of the exchanged meson, u 
and å are the usual Dirac spinor and its adjoint (Gu = 1,u = ut"), gı and ga are the coupling constants 
at the two vertices, being the Us their corresponding Dirac structures 


P,=1, Peay, Mey; eo, Lo = 8, (6) 


In the case of scalar and pseudoscalar meson exchanges, the numerator Pg of the propagator is 
just 1. For vector-meson exchange, however, it is the rank-2 tensor 
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- ut, 
e Puw = art mL (7) 


where g,, = diag(l, —1, —1, —1) is the usual Minkowski metric tensor and q, = (pi — pi), is the 
four-momentum transfer. 


In general, when all types of baryons are included, the interaction potential will be simply the sum 
of all the partial contributions 


(pıp2|V |pipa) = X (ppl Ve [P1p2) - (8) 
® 


Expanding the free Dirac spinor in terms of 1/M (where M is the mass of the relevant baryon) to 
lowest order leads to the familiar non-relativistic expressions for the baryon-baryon potentials, which 
through Fourier transformation give the configuration-space version of the interaction. The general 
expression for the local approximation of the baryon-baryon interaction in configuration space is 


vir) = Y f Cos + ma een Cis ( : d —x)rs 


mer 
P o 


+ OC, (1+ 3 : ) suo SE. (9) 
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where Co, Cos, Cis, and Cr, are numerical factors containing the baryon-baryon-meson coupling 
constants and the baryon masses, L is the total orbital angular momentum, S is the total spin, and 
Si»(f) is the usual tensor operator in configuration space, 


Si (f) = 3(01 T)(os- T) — (01:02), fF = —. (10) 


Finally, one has to remember that in the meson exchange theory all meson-baryon vertices must be 
necessarily modified by the introduction of the so-called form factors. Each vertex is multiplied by a 
form factor of the type 


A2 — m2 \™ 
F (|k = | == 11 
NEE GA (11) 
or by 
Fl) = exp ( 5E (12) 
S 242) ` 


In Eq. (11) the quantity na is usually taken equal to 1 (monopole form factor) or 2 (dipole form factor). 
The vector k denotes the three-momentum transfer, whereas A, is the so-called cut-off mass, typically 
of the order 1.2 — 2 GeV. Originally the form factors were introduced for purely mathematical reasons, 
namely to avoid divergences in the scattering equation. Nowadays our present knowledge of the (quark) 
substructure of baryons and mesons provides a physical reason for their presence. Obviously, it does 
not make sense to take the meson-exchange picture seriously in a region in which modifications due to 
the extended structure of hadrons come into play. 


2.1.2 Potential models 


As far as the potential models are concerned, they have a quite complex structure expressed via oper- 
ator invariants consistent with the symmetries of the strong interactions, i.e. translational invariance, 
Galilean invariance, rotational invariance, space-reflection invariance, time-reversal invariance, invari- 
ance under the interchange of particles 1 and 2, isospin symmetry, and hermiticity. The most widely 
known potential models are the Urbana [34] and Argonne potential [35], in which a sufficiently generic 
form of the NN potential acting between a nucleon pair 77, able to reproduce the abundance of NN 
scattering data, is expressed in operatorial form as 


18 
Vi = Kë Valrig)OG ; (13) 
u=1 


with 


A A A 


O =1, Ti" Tj, Oi: Oj, (07: 05) (Ti: Tj), Sig, Sig (Ti Tj) 
ES poses PL), 
Lite: - alm 73), (L- SY, (L- Sites - 73) - (14) 


The notation r;; = r; — rj indicates the relative position vector. The operators e; and g; are spins, 
whereas 7; and 7; are isospins (both in units of h/2), and the tensor operator is given in Eq. (10). The 
relative momentum is denoted by pi; = Di —D;; L= ri; X Pij is the total orbital angular momentum, and 
Ê? its square in the centre-of-mass system. The spin-orbit coupling enters via L-S, being $ = (0;+07)/2 
the total spin (in units of A). The first fourteen operators are charge-independent, whereas the terms 
with Gr are small and break charge independence, and they correspond to Vap(T = 1) = Van = 
Vpp, While the charge symmetry implies only that Vpn = Vp». Modern fits to very precise nucleon 
scattering data indicate the existence of charge-independence breaking, but the effect of such forces on 
the energy of nucleonic matter is much smaller than the uncertainties of many-body calculations and 
therefore can be neglected while constructing the EoS. 

Calculations show that the two-body interactions which satisfactorily reproduce NN scattering and 
deuteron properties, give the binding energies of ?H and “He systematically lower than experimental 
ones. This indicates the necessity of introducing TBFs into the nuclear Hamiltonian. Moreover, realistic 
two-body forces saturate the nuclear matter at a density significantly higher than 0.16fm >. The TBF 
can correct this, because it supplies a repulsion which increases rapidly with the growth of nucleon 
density. The TBF Maut, depends on spatial, spin and isospin coordinates of three nucleons and cannot be 
reduced to a sum of two-body interactions involving these coordinates. Widely used phenomenological 
models of TBF are dubbed as Urbana UVII and UIX [36, 37], which are expressed by the following 
decomposition 

Vige = Våg + Vor, (15) 
where the two-pion exchange part is dominant at low densities and provides the additional binding 
for the ?H and *He nuclei, whereas the IS part gives a repulsion at higher densities that is needed to 
saturate correctly the symmetric nuclear matter (SNM); this part does not contribute much in low- 
density systems such as light nuclei. 

Besides phenomenological models, microscopic approaches to the calculations of TBF are available, 
in which the same meson-exchange parameters as the two-body potentials are employed [38-42]. Al- 
though this approach has been frequently adopted in the EoS calculations of nuclear matter, the theory 
still needs to be refined, and further studies are required. 


2.2 Chiral perturbation expansion 


A different approach to the study of the baryon-baryon interactions is the one based on quark and gluon 
degrees of freedom, thus connecting the low-energy hadron physics phenomena with the underlying QCD 
structure of the baryons. However, deriving nuclear forces directly from QCD is problematic, because 
the whole hadron sector is in the non-perturbative regime, due to confinement. In fact each nucleon is, 
by itself, a complicated many-body system consisting of quarks, quark-antiquark pairs and gluons, thus 
rendering the two-nucleon problem an even more complex many-body problem. Second, the interaction 
among quarks, which is due to the exchange of gluons, is very strong at the low energies involved in 
nuclear physics processes. For this reason, it is difficult to find converging perturbative solutions. 

A new era for the theory of nuclear forces started when Steven Weinberg [20, 21] worked out an 
effective field theory (EFT) for low-energy QCD. He argued that all one needs to do is to write the 
most general Lagrangian consistent with all the properties of low-energy QCD, as this would render the 


Figure 1: Leading-order (upper diagrams) and next-to-leading-order (lower diagrams) contributions 
to the baryon-baryon interaction. Figure adapted from Ref. [52]. 


theory equivalent to low-energy QCD. A crucially important property for this discussion is the chiral 
symmetry, which is spontaneously broken, since the bare u and d quark masses in non-strange matter 
are just a few MeV. According to a theorem first proven by Goldstone, the spontaneous symmetry 
breaking implies the existence of a pseudoscalar meson, the pion. Thus, the pion plays an outstanding 
role in generating the nuclear force. Along this line, Weinberg [20, 21] proposed a scheme for including 
in the interaction a series of operators which reflect the partially broken chiral symmetry of QCD. 
The strength parameters associated to each operator are then determined by fitting the properties of 
few-body nuclear systems, i.e. deuteron and NN phase shifts. The method is then implemented in the 
framework of the Effective Field Theory (EFT) by ordering the terms according to their dependence 
on the physical parameter q/m, where m is the nucleon mass and q a generic momentum that appears 
in the Feynman diagram for the considered process. This parameter is assumed to be small and each 
term is dependent on a given power of this parameter, thus fixing its relevance. 

In this way the chiral perturbation expansion (ChPE) allows to calculate the various contributions 
to the potential systematically order by order, where each order refers to a particular power of the 
momentum. In particular, the pion-exchange term is treated explicitly and is considered the lowest- 
order (LO) term of the expansion. Furthermore, the ChPE can generate not only the force between two 
nucleons, but also many-nucleon forces in a consistent manner; in fact the TBFs so introduced are of 
higher order than the simplest two-body forces and they are treated on an equal footing. They arise first 
at next-to-next-to-leading order (N2LO) and, as a consequence, because of the hierarchy intrinsic in the 
chiral expansion, TBFs are expected to be smaller than two-body forces, at least within the range of 
validity of the expansion, whereas four-body forces appear only at next-to-next-to-next-to-leading order 
(N3LO) level, and so on. This ChPE can be used to construct NN interactions that are of reasonably 
good quality in reproducing the two-body data [43, 44]. The assumption of a small q/m parameter in 
principle restricts the applications of these forces to not too large momenta, and therefore to a not too 
large density of nuclear matter. It turns out that the safe maximum density is around the saturation 
value, ny. During the last years the NN interaction has been described to high precision using ChPE 
[45-47]. At present the nucleonic interaction has been calculated up to N4LO [48]. An exhaustive list 
of higher-order diagrams up to N3LO can be found in review papers [49, 50]. This method has been 
refined along the years and many applications can be found in the literature. 

There are very few studies of the YN interaction using ChPE in comparison to the NN case. A 
recent extension of the scheme used in Ref. [46] to the YN and the YY interactions has been performed 
by the Jiilich-Bonn-Munich group [51—54]. In the following, we briefly describe the ChPE approach to 
the baryon-baryon interaction at the two lowest orders (LO and NLO) of the chiral expansion and refer 
the interested reader to the original works of this approach for details. 


At leading order in the power counting the baryon-baryon potential consist of one pseudoscalar- 
meson exchange and of four-baryon contact terms, where each of these two contributions is con- 
strained via SU(3)-flavor symmetry (see the upper diagrams of Fig. 1). The contribution from the 
one pseudoscalar-meson exchange term is obtained from the Lagrangian density 


" m? F. 
L = (iBy"D,B — MyBB + 5 PY It Up, B}+ 3 Plu, Bl, (16) 


where the brackets denote the trace in flavor space, B is the irreducible baryon octet representation of 
SU(3)-flavor given by 

XLA + 

ye > p 
B= 3X -5+4 n (17) 
H- z — 23 


v6 
where D, is the covariant derivative, Mo is the octet baryon mass in the chiral limit, F and D are 
coupling constants satisfying the relation F + D = ga œ 1.26 with ga the axial-vector strength and 


Up = i(uld,u — 10, ul) with 
iP 
u = exp | —— |, 18 
P E ES 


being F, = 92.4 MeV the weak pion decay constant, and 


7 7 j T 
atv 7 pal 
Si = T 7] 
_K- K 29 


the SU(3)-flavor irreducible octet representation of the pseudoscalar mesons. The form of the baryon- 
baryon potentials obtained from this contribution is similar to the ones derived from the meson-exchange 
approach and in momentum space reads 


BB — (c1: q) (02 : q) 
Vope = — fri BaPSBaBaP q? + må, 


Tp, B+ Bs Bs (20) 


with fp,p»p and Topp the coupling constants of the two vertices, Mps the mass of the exchanged 
pseudoscalar meson, q the transferred momentum, and Tz, p, , p,p, the corresponding isospin factor. 

The contribution from the four-baryon contact interactions can be derived from the following mini- 
mal set of Lagrangian densities 


L=C (BB TB ys )a) , (21 
L> = C? (Ba(T:B)aBi(T:B)o) D (22 
£ =CHBA(T¡B)a)(Br(T¡B)o) (23 


Here, the labels a and b are the Dirac indices of the particles and I; denotes the five elements of the 
Clifford algebra, Tı = 1,19 = yt, I'3 = o"",T4 = yy’, T's = y, which are actually diagonal 3 x 3 
matrices in the flavor space. In LO these Lagrangian densities give rise to six independent low-energy 
coefficients (LECs): Ch, C1, C2, C2,C2 and C2, where the labels C and S refer to the central and 
spin-spin parts of the potential, respectively. The LO contact potentials for the different baryon-baryon 
interactions resulting from these Lagrangians have the general form 


Vio = CE" + 03% a1 - 02, (24) 
where the coefficients CEP and CH? are linear combinations of CL, C1, C2, C2, C3 and C3. 
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At next-to-leading order (see the lower diagrams of Fig. 1), the contact terms read 


Vara = Cig? + Cok? + (Caq? + Cık?)oı -O2 + Oso, + 02) : (q x k) 
t Cela : 01)(g 02) + C7(k - o1)(k 07) + Oslo, - 02) (q x k), (25) 


where C; (i = 1,--- ,8) are additional LECs. The momenta q and k are defined in terms of the initial p 
and final p' baryon momenta in the center-of-mass frame as q = p' — p and k = (p+ p')/2, respectively. 
The expressions for the two-pseudoscalar meson-exchange contributions are rather cumbersome and we 
refer the interested reader to the original works for details (see e.g. Ref. [52]). 

The baryon-baryon potentials constructed in this way are then inserted in the Lippmann-Schwinger 
equation which is regularized with a cut-off regulator function of the type 


p + p^ 
Al 


Flat) = exp ( (26) 
in order to remove high-energy components of the baryon and pseudoscalar meson fields. The cut-off 
A is usually taken in the range 450 — 700 MeV. 


2.3 Resonating-group method 


The resonating-group method (RGM) has been used to develop a further approach to the study of the 
NN interaction [55-58]. In this approach, inspired by the QCD theory of strong interaction, the quark 
degree of freedom is explicitly introduced and the NN interaction is constructed from gluon and meson 
exchange between quarks, the latters being confined inside the nucleons. Due to the RGM formalism, 
the resulting interaction is highly non-local and contains a natural cut-off in momentum. The most 
recent model, named fss2 [55, 59], reproduces closely the experimental phase shifts, and fairly well the 
data on the few-body systems, e.g. the triton binding energy is reproduced within 300 keV. This model 
has been used in nuclear matter calculations performed with the non-relativistic BHF approach, and 
in particular it has been found that it is able to reproduce correctly the nuclear matter saturation 
point without the TBF contribution [60]. This finding is quite relevant, and surely deserves further 
investigation. 


2.4 Renormalization-group methods 


A further class of NN interactions is based on renormalization group (RG) methods (see e.g. [13, 61] for 
a complete review). The main effect of the hard core in the NN interaction is to produce scattering to 
high momenta of the interacting particles. A possible way to soften the hard core from the beginning is 
by integrating out all the momenta larger than a certain cut-off A and “renormalize” the interaction to 
an effective interaction Vowk in such a way that it is equivalent to the original interaction for momenta 
q < A. This results in a modified Lippmann-Schwinger equation with a cutoff-dependent effective 
potential View i 


2 f^ owk(k’, QT (q, k; k? 
T(k',k; k’) = Viowklk’, k) + =P | ie ek q) (q ) . (27) 


oi 
0 

By demanding dT(K', k; k?)/dA = 0, an exact Renormalization Group flow equation for Vowk can be 
obtained 


dViowk(Kk!, k) — 2 Monk, A)T(A, k; A?) 

dA A 1 — k?/A? 
Integrating this flow equation one can obtain a phase-shift, energy independent, soft (i.e. without hard 
core) and hermitian low-momentum potential Vowk. The Viowx interaction turns out to be much softer, 


(28) 
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since no high momentum components are present and, as a consequence, three- and many-body forces 
emerge automatically from a pure two-body force. The short range repulsion is replaced by the non-local 
structure of the interaction. The cut-off A is taken above 300 MeV in the laboratory, corresponding to 
relative momentum q = 2.1fm”!, that is the largest energy where the experimental data are established. 
The fact that Viowx is soft has the advantage to be much more manageable than a hard-core interaction, 
in particular it can be used in perturbation expansion and in nuclear structure calculations in a more 
efficient way [61, 62]. 


Following the same idea that in the NN case made possible to calculate a “universal” effective 
low-momentum potential Vowx by using renormalization-group techniques, in [63] this method was 
generalized to the YN sector. Unfortunately, contrary to the NN case there exist only few YN scattering 
data and hence the YN interaction is not well constrained. It was found (see Figs. 1-6 of Ref. [63]) 
that the YN phase shifts have approximately the same shape but different heights, and the diagonal 
matrix elements differ for lower momenta, although they collapse for momenta near the cut-off. In 
conclusion, however, one can still say that in general the results seem to indicate a similar convergence 
to an “universal” softer low-momentum YN interaction as for the NN case. 


2.5 Baryon-baryon interactions from lattice QCD 


Recently, a further possibility of constructing the baryon-baryon interaction based on lattice QCD has 
been explored, see [64-66] for a review. Lattice QCD, from which one should be able in principle to 
calculate the hadron properties directly from the QCD Lagrangian, is however extremely expensive from 
the numerical point of view and current simulations can be performed only with large quark masses [07]. 
In fact, an accurate simulation has to be made on a fine grid spacing and large volumes, thus requiring 
high-performance computers. Nonetheless, a big progress to derive baryon-baryon interactions from 
lattice QCD has been made by the NPLQCD [68] and HALQCD [69] collaborations. It is worthwhile to 
point out, however, that there exist some discrepancies between these two collaborations regarding the 
methods employed in their studies. The NPLQCD collaboration combines calculations of correlation 
functions at several light-quark-mass values with the low-energy effective field theory. This approach 
is particularly interesting since it allows to match lattice QCD results with low-energy effective field 
theories providing the means for first predictions in the physical quark mass limit. The HALQCD 
collaboration, on the other hand, follows a method to extract the different baryon-baryon potentials 
from the Nambu-Bethe-Salpeter wave function measured on the lattice. Recently this collaboration 
managed to approach the region of physical masses, obtaining results for various NN, NY and YY 
interaction channels at a single value of the lattice volume and of the lattice spacing [70-72]. 


Using the NN forces thus obtained, the properties of nuclei such as *He, 16O and Ca and the 
EoS of nuclear matter were investigated [73, 74], finding that these nuclei and the symmetric nuclear 
matter are bound at a quark mass corresponding to a pseudo-scalar meson (pion) mass of 469 MeV. The 
obtained binding energy per nucleon has a uniform mass-number A dependence, which is qualitatively 
consistent with the Bethe-Weizsácker mass formula, thus demonstrating that the HALQCD strategy 
works well to investigate various properties of atomic nuclei and nuclear matter starting from QCD. 
We would also like to note that in the strangeness sector the NPLQCD collaboration has been able to 
determine the binding energies of light hypernuclei including å He, 4 He and å, He [75]; to compute the 
magnetic moment of the baryon octet [76] and to constrain the interactions of two-baryon octets at the 
SU(3)-flavor symmetric-point [77]. 
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3 Theoretical approaches to the EoS of the core 


Once the interaction between nucleons is established, one can try to solve the many-body problem for 
nuclear matter. Since there are currently no ab initio QCD calculations of dense matter available at 
the thermodynamic conditions typical of compact stars, one has to rely on interaction models, which 
can be constrained by laboratory measurements and NS observations. In the following we discuss 
mostly theoretical approaches based on nucleonic degrees of freedom, whereas later on we will treat 
approaches where hyperons or even quarks are included, since they are essential at large densities. The 
currently used many-body methods are divided into two classes, i.e. i) Ab-initio microscopic methods 
which start from bare two- and three-nucleon interactions able to reproduce nucleon scattering data and 
properties of bound few-nucleon systems. The main drawback of the ab-initio methods is that they can 
be employed only for the description of homogeneous matter in the NS core, and not for the clustered 
matter typical of the NS crust; ii) Phenomenological approaches which use effective interactions with a 
simple structure dependent on a limited number of parameters, which are fitted to different properties 
of finite nuclei and nuclear matter. The main problem of these approaches is their extrapolation to 
high-density conditions, because those models are fitted to the ground-state properties of finite nuclei. 

In the following we will briefly discuss all of them, referring the interested reader to Refs. [78, 79] 
and the quoted references in each subsection. 


3.1 Ab-initio microscopic methods 
3.1.1 (Dirac) Brueckner-Hartree-Fock 


The Bethe-Brueckner-Goldstone (BBG) many-body theory is based on the re-summation of the per- 
turbation expansion of the ground-state energy of nuclear matter [80, 81]. The original bare NN 
interaction is systematically replaced by the so-called G-matrix, an effective interaction that describes 
the in-medium scattering processes, and that takes into account the effect of the Pauli principle on the 
scattered particles. A further essential ingredient is the in-medium single-particle potential U(k) felt by 
each nucleon, where k is the momentum. The key point is the solution of the Bethe-Goldstone integral 
equation for the G-matrix, which can be written as 


[1 — Or(k3)] [1 — Or (k4)] 


kk IG kek 29 
W — ek, + eg, + (k3k4|G(w)|k3ka) , (29) 


(kik2|G(w)|k3ka) = (kıka|V |kaka) +) (kıka|V|k3K4) 


Kb kl 


where V is the bare NN interaction, w is the starting energy, the two factors [1 — Or(k)] force the 
intermediate momenta to be above the Fermi momentum (“particle states”), the single-particle energy 
being ep = k?/2m+U(%), with m the particle mass, and the summation includes spin-isospin variables. 
The main feature of the G-matrix is that it does not have the hard core of the original bare NN 
interaction and it is defined even for bare interactions with an infinite hard core. The introduction 
and choice of the in-medium single-particle potential are essential to make the re-summed expansion 
convergent, and is calculated self-consistently with the G-matrix itself in order to incorporate as much 
higher-order correlations as possible. The resulting nuclear EoS can be calculated with good accuracy at 
the so-called Brueckner-Hartree-Fock (BHF) level (“two hole-line” approximation) with the continuous 
choice for the single-particle potential, the results in this scheme being quite close to the calculations 
which include also the three hole-line contribution [82-84]. 

It should be noted that all non-relativistic many-body approaches fail to reproduce the correct 
saturation point of nuclear matter. For this purpose, TBFs are introduced in the calculations, and are 
reduced to a density dependent two-body force by averaging over the generalized coordinates (position, 
spin, and isospin) of the third particle. In the BHF calculations for nuclear matter, the so-called Urbana 
model is often introduced as TBF, and this produces a shift in the binding energy of about +1 MeV 
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and of —0.01fm ? in density. The problem of such a procedure is that the TBF is dependent on 
the two-body force. The connection between two-body and TBFs within the meson-nucleon theory of 
nuclear interaction is extensively discussed and developed in [39, 40, 42]. At present the theoretical 
status of microscopically derived TBFs is still quite rudimentary; however, a tentative approach has 
been proposed using the same meson-exchange parameters as the underlying NN potential. Results 
have been obtained with the Argonne Vis [35], the Bonn B [85], and the Nijmegen 93 potentials [41, 42]. 
However, the role of TBF also depends on the model adopted for TBF. In fact, if the NN potential 
is based on a realistic constituent quark model [60] then the role of TBF is strongly reduced. On the 
other hand, if we consider a new class of chiral-inspired TBF, then the BHF calculations are not able 
to reproduce simultaneously the correct saturation point and the properties of three- and four-nucleon 
systems [86]. 

The relativistic BHF formalism, i.e. the Dirac-Brueckner approach [87], has been developed in 
analogy with the non-relativistic case, where the two-body correlations are described by introducing 
the in-medium relativistic G-matrix. The novel and most striking feature of the DBHF theory is its 
ability to describe the saturation properties of nuclear matter, a fundamental aspect which reflects 
the saturating nature of the nuclear force. The DBHF method contains important relativistic features 
through the description of the nuclear mean field in terms of a scalar and a vector component, strong 
and of opposite sign. In their combination, they provide an explanation for the binding of nucleons 
and the spin-orbit splitting in nuclear states. The main relativistic effect is due to the use of the 
spinor formalism, which has been shown [88] to be equivalent to introducing a particular TBF, the 
so-called Z-diagram, which is repulsive and consequently gives a better saturation point than the BHF 
method, and actually including in BHF only these particular TBFs, one gets results close to DBHF 
calculations [89]. In conclusion, the EoS calculated within the DBHF method turns out to be stiffer 
above saturation than the ones calculated from the BHF+TBF method. A further point of difference 
is that the relativistic DBHF produces a superluminal EoS at higher densities than the BHF approach. 
The reader is referred to Ref. [90] for a relatively recent review of the DBHF method and a variety of 
applications to both nuclear matter and nuclei. 


3.1.2 Approaches based on the variational method 


The variational approach to the many-body problem is based on the Ritz-Raleigh variational principle, 
according to which the trial ground state energy 


(vH) 
Errial — (UD) , (30) 


calculated from the system's Hamiltonian H with a trial many-body wave function V, gives an upper 


bound for the true ground state energy of the system. In the variational method one assumes that the 
ground state wave function V can be written as 


V(ri 12, ......) = [TF (rig) (ri, ra...) , (31) 


where ® is the unperturbed ground-state wave function, properly antisymmetrised, and the product runs 
over all possible distinct pairs of particles. The correlation factors f, which are intended to transform 
the uncorrelated wave function ® to the correlated one, are determined by the Ritz-Raleigh variational 
principle, i.e. by assuming that the mean value of the Hamiltonian reaches a minimum: 


5 (ww) _ 
5f (UD) 0. (32) 
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Therefore the main task in the variational method is to find a suitable ansatz for the correlation 
factors f. Several different methods exist for the calculation of f, e.g. in the nuclear context the 
Fermi-Hyper-Netted-Chain (FHNC) [91, 92] calculations have been proved to be efficient. 

For nuclear matter it is necessary to introduce channel-dependent correlation factors, which is equiv- 
alent to assume that f is actually a two-body operator 129 One then assumes that P can be expanded 
in the same spin-isospin, spin-orbit, and tensor operators appearing in the NN interaction [93, 94]. Due 
to the formal structure of the Argonne NN forces, most variational calculations have been performed 
with this class of NN interactions, often supplemented by the Urbana TBFs. The best known and most 
used variational nuclear matter EoS is the one of Akmal-Pandharipande-Ravenhall (APR) [95]. Many 
excellent review papers exist in the literature on the variational method and its extensive use for the 
determination of the nuclear matter EoS, e.g. [92]. 

Among other methods based on the variational principle and widely used in nuclear physics, we men- 
tion the coupled-cluster theory, proposed in [96, 97], in which the correlation operator is represented in 
terms of the cluster operators. In recent nuclear matter calculations with chiral NN interactions [98, 99] 
and also in nuclear structure calculations [100, 101], the method has been proven to be successful. The 
variational Monte Carlo (VMC) approach is also commonly used, especially for light nuclei, including 
two and three-body correlations [102], but the EoS of homogeneous nuclear matter is hard to obtain, 
due to the large number of nucleons which increases the computational effort [94, 103]. 


3.1.3 Self-consistent Green's function 


Another way to approach the many-body problem is through the many-body Green's functions for- 
malism [104]. In this approach one performs a diagrammatic analysis of the many-body propagators 
in terms of free one-body Green's functions and two-body interactions. The perturbative expansion 
results in an infinite series of diagrams, among which one has to choose those which are relevant for the 
considered physical problem. Depending on the approximation, one can either choose a given number 
of diagrams or sum an infinite series of them, in analogy with the BBG approach. In the description 
of nuclear matter, the method is conventionally applied at the ladder-approximation level, which en- 
compasses at once particle-particle and hole-hole propagation, and this represents the main difference 
with respect to the G-matrix, where only particle-particle propagators are included. At a formal level, 
the comparison between the BBG and the SCGF approaches is not straightforward. Even though both 
approaches arise from a diagrammatic expansion, the infinite subsets of diagrams considered in the two 
approaches are not the same, and the summation procedures are also somewhat different. Whereas the 
BHF formalism in the continuous choice can be derived from the ladder SCGF formalism after a series 
of approximations, this is not the case for the full BBG expansion. In principle, if both BBG and SCGF 
were carried out to all orders, they should yield identical results. BBG theory, however, is an expansion 
in powers of density (or hole-lines), and the three-hole line results seem to indicate that it converges 
quickly. The error in the SCGF expansion is more difficult to quantify, as one cannot directly compute 
(or even estimate) which diagrams have to be included in the expansion. 

The energy per particle of nuclear matter is obtained within the SCGF approach through the 
Galitskii-Migdal-Koltun sum-rule [105, 106]: 


STENE naro, o 


where y = 4(2) is the spin-isospin degeneracy of nuclear (neutron) matter, p is the total density and 
f(w) =[1+e2/T]7 is the Fermi-Dirac distribution. The key quantity of the approach is the one-body 
spectral function A(k,w) which, in few words, represents the probability density of removing from or 
adding to the system a nucleon with momentum k and energy w. The single-particle spectral function 
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gives access to the calculation of all the one-body properties of the system and it can be obtained from 
the imaginary part of the one-body propagator G(k, w), 


A(k,w) = git (34) 


or, alternatively, from the proper or irreducible self-energy 2(k,w) as 


Im X(k,w) 


459) =~ — E -Rez(k o) n Sto 


(35) 


The computational implementation of the SCGF method requires: (1) to calculate the effective 
interaction (T-matrix) describing the in-medium scattering of two nucleons, (2) to extract the self 
energy X(k,w) and (3) to obtain the one-body propagator G(k,w) by solving, with the self-energy, the 
Dyson equation, which is then inserted in the scattering equation, repeating these three steps until a 
self-consistent solution is achieved. 

Also for the SCGF method the inclusion of TBFs is essential. So far TBFs were not included in 
the ladder approximation, however a method has been developed recently in [107], and applied to SNM 
using chiral nuclear interactions. TBFs are included via effective one-body and two-body interactions, 
and are found to improve substantially the saturation point [105]. 

One has to note that because of the well-known Cooper instability [109], through which a fermionic 
many-body system with an attractive interaction tends to form pairs at the Fermi surface, low-temperature 
nuclear matter is unstable with respect to the formation of a superfluid or superconducting state. The 
Cooper instability shows up as a pole in the 7-matrix when the temperature falls below the critical 
temperature T, for the transition to the superfluid/superconducting state. Therefore current calcula- 
tions are often performed at temperatures above T, and extrapolated to zero temperature, see e.g. [110]. 
Further details on the SCGF method and on its applications to nuclear problems can be found, e.g. in 
[111-113]. 


3.1.4 Chiral effective field theory (xEFT) approach 


High-precision nuclear potentials based on the ChPE [45] are nowadays very popular, since ChPE allows 
to link nuclear physics with QCD. In particular, the treatment of many-nucleon forces is very important, 
those being mostly relevant in nuclear matter. In the last years some effort has been devoted to the 
extension of chiral perturbation theory to nuclear matter calculations, i.e. developing an effective field 
theory (EFT) for nuclear matter. In this case another scale appears, kr/M, being M the nucleon mass 
and kr the Fermi momentum, approximately given as kp = 263 MeV at saturation which is smaller than 
a typical hadron scale. In the chiral limit it is then natural to expand in kr/M, and this expansion can 
be obtained from the vacuum CHPE [114]. For nuclear matter the correction thus obtained with respect 
to the vacuum diagrams gives a direct contribution to the EoS of nuclear matter, and this correction 
is clearly proportional to a power of kr /M. In this case a cut-off must be introduced, and its tuning 
allows to obtain a saturation point and compressibility in fair agreement with phenomenology. Along 
the same lines more sophisticated approaches can be developed, where the many-nucleon interactions 
built in vacuum are directly used in nuclear matter calculations. In this case the ChPE is used in 
conjunction with the EFT scheme (xEFT). 

In recent years, xEFT has been used for studying nuclear matter within various theoretical frame- 
works like many-body perturbation theory [115-118], SCGF framework [108], in-medium chiral per- 
turbation theory [119], the BHF approach [120, 121] and quantum Monte Carlo methods [122-124]. 
Several reliable calculations have been performed up to twice the saturation density ny, beyond which 
uncertainties were estimated by analysing the order-by-order convergence in the chiral expansion and 
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the many-body perturbation theory [125]. Variations in the resolution scale [126] and low-energy con- 
stants appearing in the two-nucleon and three-nucleon forces were systematically explored [127], finding 
that the theoretical uncertainty band grows rapidly with the density beyond ny, due to the missing 
third-order terms at low densities and higher-order contributions in the chiral expansion. This has 
consequences not only for the EoS, but also for the symmetry energy at saturation density, So, and the 
slope parameter L, as discussed in [128]. 


3.1.5 Quantum Monte Carlo methods 


Quantum Monte Carlo (QMC) methods are very successful in solving the nuclear many-body problem 
non-perturbatively and with controlled approximations, thus making QMC methods quasi-exact in 
describing the ground state of fermionic systems. Among the most successfully described, we mention 
the liquid ?He, but also bosons, like atomic liquid “He. Moreover they have been applied in studies 
of nuclear matter and light nuclei [94, 129]. Several implementations of QMC methods have been 
developed over the years, e.g. Green’s function Monte Carlo (GFMC) [130], or auxiliary-field diffusion 
Monte Carlo (AFDMC) [131], which differ in the treatment of the spin and isospin degrees of freedom. 
In particular, the AFDMC method has been used extensively to study nuclear matter for astrophysical 
applications [124, 131, 132]. 

The main idea of Quantum Monte Carlo methods is to stochastically solve the many-body Schrodinger 
equation to extract the ground state of a system, by evolving a given trial wave function of the many- 
body system, Vy, in imaginary time 7 = it. Monte Carlo sampling is then used to evaluate all possible 
configurations until convergence is reached. However, the accuracy of the different QMC versions is lim- 
ited by the so-called fermion sign problem [133], for which different approximations are adopted [134]. 
In fact, for fermionic systems the wave function is antisymmetric and contains several changes in sign. 
Hence, the integrands in the QMC integrals are highly oscillatory thus producing very large statistical 
uncertainties, so that no information can be obtained from the calculation. One possible way to cure 
the problem consists in splitting the wave function space into regions of positive and negative wave 
functions, defining a nodal surface at which the wave function changes sign. Generally, configurations 
that cross the nodal surface are removed from the evolution. This approximation, called fixed-node 
approximation [133], has been generalized to the constrained-path method [135, 136] for complex wave 
functions. 

When calculating nuclear matter, one typically simulates N particles in a cubic box with size L, 
where L is determined in such a way that the number density n in the box reflects a chosen value, 
L=(N/n)">. To probe the thermodynamic limit, the particle number N has to be chosen sufficiently 
large. However, due to growing computational costs associated with large particle numbers, for neutron 
matter one typically chooses N = 66 (33 spin up and 33 spin down neutrons), thus finding results 
close to the thermodynamic limit. Those approximations seriously limit the potentiality of the QMC 
methods. Moreover, we notice that in spite of its recent progress, it is not yet possible to perform 
GFMC and AFDMC calculations with the Argonne Vig potential, mainly due to technical problems 
associated with the spin-orbit structure of the interaction and the trial wave function, which induce very 
large statistical errors. In order to overcome this problem, the full operatorial structure of current high- 
quality NN potentials has been simplified and more manageable NN potentials have been developed 
containing less operators with readjusted parameters. In particular, we mention the Vi, Vg and Vj 
potentials [137, 138], eventually supplemented with the Urbana TBFs. Recently, a local chiral potential 
has been developed [122], which is well suited for QMC techniques. 

Modern computer technology has allowed the extension of the QMC method to nuclear systems, 
which have more complicated interactions and correlation structures. It has to be noticed that the com- 
puting time increases exponentially with the number of particles, which limits the number of nucleons 
considered by GFMC up to 16 neutrons. The largest nucleus considered is !?C. A recent comparison 
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has demonstrated that both methods give very close results for neutron drops with N < 16 [139]. 


3.2 Phenomenological approaches 
3.2.1 Non-relativistic energy density functionals 


Phenomenological approaches make use of effective interactions instead of bare ones to treat dense 
matter, and mostly rely on the energy density functional (EDF) theory, which is able to recast the 
complex many-body problem of interacting particles into an effective independent particle approach 
[140]. In the EDF theory, the total energy of the system is usually expressed as a functional of the 
nucleon number densities, the kinetic energy densities, and the spin-current densities, but the exact 
form of the functional itself is not known a priori. Therefore, one has to rely on phenomenological 
functionals, which depend on a certain number of parameters fitted to reproduce some properties of 
known nuclei and nuclear matter, as well as ab-initio calculations of infinite nuclear matter. A different 
approach to construct a phenomenological EoS is to use a purely parametrized EoS, without considering 
the main properties of the NN interaction. An example is given by a metamodel for the nucleonic EoS, 
which consists in a Taylor expansion around the saturation density of SNM in terms of the empirical 
parameters [141], and has been applied to study NS global properties [142]. 

Non-relativistic approaches usually start from an Hamiltonian Å — T +V for the many-body sys- 
tem, where T = >>, Di däm: is the kinetic term (f being the momentum operator and m; the mass of 
the species 7) and V is the potential term. The latter accounts for the two-body (pseudo)potential, 
that allows one to incorporate physical properties like effective masses, but TBFs can also be included 
explicitely [143]. A very popular scheme is based on the Skyrme-type effective interactions, which are 
zero-range density-dependent interactions and allow for fast numerical computations. Since the pioneer 
work of Skyrme [144], several extensions have been proposed, which accurately reproduce experimen- 
tally measured properties of finite nuclei [145-149] and are applied to NSs. We mention that many 
Skyrme models have been tested against several nuclear-matter constraints in Ref. [150], but most of 
the constraints are known with large error bars, particularly the symmetry energy coefficients, and 
therefore it is not possible to rule out some models on this basis. Besides Skyrme models, finite-range 
(density-dependent) interactions have been derived from the Gogny interaction [151], which are less 
widely used in astrophysics because of the larger numerical efforts [152, 153]. 

In addition to the Skyrme and Gogny effective interactions, new EDFs have been recently constructed 
on the basis of the Kohn-Sham density functional theory [154, 155]. These Barcelona-Catania-Paris(- 
Madrid) (BCP and BCPM) EDFs have been derived by introducing in the functional results from 
microscopic nuclear and neutron-matter BHF calculations, thus yielding a very good description of 
properties of finite nuclei and NS masses with a reduced number of parameters. 


3.2.2 Relativistic mean-field (RMF) models 


Relativistic mean-field (RMF) models have been constructed on the basis of the quantum hadrodynamics 
(QHD) framework [156], a field-theoretical formalism where nucleons are represented by four-component 
Dirac spinors, and the NN interaction is modeled by exchange of mesons. A nucleus is thus described 
as a system of Dirac nucleons whose motion is governed by the Dirac equation. RMF models have been 
successfully employed in nuclear structure, to describe both nuclei close to the valley of stability and 
exotic nuclei [157]. 

The common starting point in RMF models is an effective Lagrangian L = Lyye + Lmes + Lint, 
where the different terms account for the nucleon, the free mesons and the interaction contribution, 
respectively. The isoscalar scalar o-meson and the isoscalar vector w-meson mediate the long-range 
attraction and short-range repulsion of the nuclear interaction, respectively, in SNM. Isovector mesons 
(like the isovector vector p-meson and the isovector scalar ó-meson) need to be included as well to 
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treat neutron-proton asymmetric systems. The interaction term depends on the nucleon-meson cou- 
pling constants that are usually determined by fitting nuclei or nuclear-matter properties. From the 
Lagrangian above, the field equations for nucleons and mesons are derived, and they are solved self- 
consistently, usually in the relativistic mean-field approximation. However, the correct description of 
nuclear matter and finite nuclei requires the extension of that simple Lagrangian in order to include a 
medium-dependent effective interaction. This effect can either be introduced by including non-linear 
(NL) meson self-interaction terms in the Lagrangian, or by assuming an explicit density dependence 
(DD) for the meson-nucleon couplings, which gives rise to the so-called rearrangement contributions 
mostly important for the thermodynamic consistency of the model. 

The former approach has been employed in constructing several phenomenological RMF interactions, 
like the widely used NL3 [158], PK1, PKIR [159] and FSUGold [160]. In the second approach, the 
functional form of the density dependence of the coupling constants can be derived by comparing results 
with the microscopic DBHF calculations or it can be fully phenomenological, with parameters adjusted 
to experimental data (DD-RMF models, see Refs. [159, 161, 162]). We notice that all those models do 
not explicitly take into account the antisymmetrization of the many-body wave function. This has been 
included in relativistic Hartree-Fock (RHF) or Hartree-Fock-Bogoliubov models accounting for pairing, 
see Refs. [163-165]. 

Due to the mean-field character of the model and the charge independence of the strong interaction, 
the extension of the model to the full baryon octet is slightly more difficult than the simple nuclear- 
matter case. In Ref. [166], an ensemble of 263 RMF models was tested against nuclear matter properties 
up to about 3 times the saturation density, also making use of heavy-ion collision data. Only a limited 
set was found to be consistent with the constraints taken into account. This same set has been further 
used in Ref. [167], where the Love number and the corresponding tidal deformabilities show very good 
agreement with the data from the GW170817 merger event. 

A further RMF model is the quark-meson-coupling model, which incorporates the internal quark 
structure of baryons, where nucleons are treated as bound states of three quarks and interact via meson 
exchange, pions included. This model has been applied to study NS properties, e.g. in [168]. 

Both Skyrme and RMF based approaches turn out to be very popular in astrophysical applications, 
i.e. core-collapse supernovae and binary NS mergers simulations, where a wide range of densities, tem- 
peratures and charge fractions, describing both clustered and homogeneous matter, is needed. This 
range is covered by the so-called “general-purpose” EoSs, which are described in Sec. 3.4. However, 
we like to point out that phenomenological models have to be considered with some care, as discussed 
in Ref. [169], where a comparison between some phenomenological EoSs with those derived within the 
chiral EFT interactions up to next-to-next-to-next-to-leading order (N3LO) shows some inconsistency 
with the chiral pure neutron matter (PNM) band, although this comparison is limited to a density 
range up to p = 0.165 fm ?. This point is further illustrated in Sec. 3.3. 


3.3 Comparing ab-initio and phenomenological approaches 


We now turn to discuss the main differences in the results obtained for SNM and PNM, employing 
microscopic and phenomenological EoS models. The main results are summarized in Fig. 2. In the 
upper (central) panels we display the binding energy per particle as a function of the nucleon density 
obtained for SNM (PNM) with a set of microscopic and phenomenological approaches and NN potentials, 
plotted respectively in the left and right columns. In the lower panels we display the symmetry energy 
for the same EoSs. 

As far as microscopic approaches are concerned, we adopt several BHF EoSs based on different 
NN potentials, namely the Bonn B (BOB) [87], the Nijmegen 93 (N93) [26, 170], and the Argonne 
Vis (V18) [35]. In all those cases, the two-body forces are supplemented by nucleonic TBFs, either 
phenomenological or microscopic [38, 39, 41, 42, 171]. In the phenomenological case the corresponding 
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Figure 2: Binding energy per nucleon of symmetric matter (upper panels), pure neutron matter (central 
panels) and the symmetry energy (lower panels) obtained with different EoSs. See text for details. 


EoS employs the Urbana model for TBF and is labelled UIX. Within the same theoretical framework, 
we also studied an EoS based on a potential model which includes explicitly the quark-gluon degrees 
of freedom, i.e. FSS2, which reproduces correctly the saturation point of SNM and the binding energy 
of few-nucleon system without the need of introducing TBF. In the following we use two different EoS 
versions labelled respectively as FSS2CC and FSS2GC. Moreover, we compare these BHF EoSs with 
the often-used results of the DBHF EoS [172], which employs the Bonn A potential, the APR EoS 
[95] based on the variational method and the Argonne Via potential, and a parametrization of a recent 
Auxiliary Field Diffusion Monte Carlo (AFDMC) calculation [173]. 

As far as phenomenological approaches are concerned, we show results obtained for a sample of seven 
Skyrme forces, namely, GS and Rs [174], SLy4 [175] of the Lyon group, the old SV [176], SkI4 [177] of the 
SkI family, SkMP [178], and SkO [179]. We also consider two types of RMF models, and in particular 
we adopt models with constant meson-baryon couplings described by the Lagrangian density of the 
nonlinear Walecka model (NLWM), and models with density-dependent couplings, hereafter referred 
to as density-dependent models (DDM). Within the first type, we consider the models GM1 and GM3 
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Figure 3: The binding energy of PNM (upper panels) for microscopic EoSs (left panels) and phe- 
nomenological ones (right panels). The light-orange band represents the predictions of xEFT theory 
reported in Ref. [134]. In the lower panels the symmetry energy is displayed, along with currently avail- 
able limits from heavy ion collisions (green, blue and gray bands), and the IAS calculations (orange 
contour). See text for details. 


[180]. For the DDM, we consider the models DDME1 and DDME2 [181], and TW99 [182]. A further 
phenomenological RMF EoS, the SFHo EoS [183], has been used for comparison. The latter belongs to 
the class of the so-called general-purpose EoSs, and will be illustrated in Sec. 3.4. 


Turning back to the discussion of Fig. 2, we see that in both PNM and SNM all approaches yield sim- 
ilar results up to about twice the saturation density, and diverge at larger density. For the microscopic 
approaches of BHF type, this is certainly due to the different high-density behaviour of TBFs, either 
microscopic or phenomenological. Also three-body correlations make a difference. For phenomenolog- 
ical models, a similar spread at high density can be observed. Indeed, these models are characterized 
by parameters which are fitted on experimental data approximately known around saturation density, 
therefore their behaviour at high density, where no experimental data are avilable, can produce large 
differences. In the lower panels of Fig. 2, we display the symmetry energy as a function of the density. 
We observe a monotonically increasing behaviour for all the considered EoSs, and strong divergencies 
at large density, which depend on either the nucleonic interaction or the many-body scheme; those are 
important for determining structure and composition of NSs, and also their cooling mechanism. 


It is worthwhile to perform a comparison of the above discussed EoSs with the currently available 
constraints. In Fig. 3 (upper panels) we compare with the results obtained for PNM within the chiral 
EFT interactions up to N3LO order [169, 184], and shown by the orange-shaded band. Please notice that 
this comparison is limited to a density range up to p = 0.34 fm ?. Whereas the microscopic calculations 
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fall well inside the chiral band, some inconsistency is evident for the phenomenological EoSs, as already 
found in Ref. [169]. In particular, we find that the SLy4, BSk22, BSk24, BSk26 and SFHo are in 
agreement with the chiral calculations; TW99, DDME1 and DDMEZ2 are not at all compatible, and 
the remaining EoSs are only partially compatible with the chiral PNM results. We stress that neutron 
matter is an ideal laboratory for nuclear interactions derived from chiral effective field theory since all 
contributions are predicted up to N3LO in the chiral expansion, including three-body forces [118]. By 
comparing with results obtained within the self-consistent Green’s function theory (SCGF) using the 
same three-body forces, predictions for the EoS of neutron matter at zero temperature at density slightly 
above the saturation density have been improved, including estimates of theoretical uncertainties for 
astrophysical applications. A recently proposed Bayesian machine-learning method allowed to improve 
results up to twice nuclear saturation density for the binding energy and pressure of neutron matter, as 
well as for the nuclear symmetry energy and its slope, which turn out to be consistent with experimental 
constraints [184] at saturation density. 

In the lower panels of Fig. 3, we display the symmetry energy over the same density range. Results 
are compared with available experimental data, which are displayed by the shaded areas. In particular, 
the grey area represents the diffusion data of HICs, the green area includes the data obtained by the 
FOPI-LAND collaboration [185] on the collective flow, and the blue area is the experimental region 
checked by the ASY-EOS collaboration [186]. The full orange contour shows the results on the isobaric 
analog states (IAS), obtained in Ref. [187]. We see that most of the considered EOSs are compatible 
with experimental data up to around saturation density, whereas for larger densities some EOSs tend 
to predict smaller values of the symmetry energy, below the experimental areas. This is a clear sign 
of discrepancy, which results in a much larger difference at larger values of the baryon density, as the 
ones characterizing the inner core of a NS. We notice that the inferred constraints are model dependent, 
since the data interpretation requires theoretical simulations. Those discrepancies give rise to different 
predictions for NS structure and cooling behaviour, as will be shown later. 


3.4 Finite-temperature EoSs 


The EoS at finite temperature of asymmetric nuclear matter plays a major role in the latest stage of the 
SN collapse [12, 13, 188] and binary NS mergers [16, 189], because it determines the final evolution of a 
transitory state to either collapse to a black hole or to the formation of a NS. A few finite-temperature 
nuclear EoSs for astrophysical simulations are now available [182, 183, 190-195], and the predictions for 
the effects of temperature on stellar stability are sometimes conflicting: RMF models usually predict 
increasing stability (maximum mass) with temperature, whereas BHF results indicate a slight reduction 
of the maximum mass. There is therefore a possible mass range of the hot rotating transitory metastable 
object that depends on the finite-temperature EoS. 

As far as microscopic calculations of the nuclear EoS at finite temperature are concerned, we mention 
in particular the non-relativistic Brueckner-like calculations [80], where the formalism by Bloch and De 
Dominicis [196] was followed. Those calculations confirmed that the resulting EoS for SNM shows a 
typical Van der Waals behaviour, which entails a liquid-gas phase transition, with a definite critical 
temperature found to be around T, = 18 — 20 MeV, very similar to the Friedman and Pandharipande 
findings [197] based on the variational method. The main drawback in this approach is that the 
Hugenholtz-Van Hove (HVH) theorem is not satisfied. In other words, the pressure p calculated from 
the thermodynamic relation p = —f + up, being f the free energy density, u the chemical potential 
and p the number density, does not coincide with the pressure calculated from p = —Q/V, Q being 
the grand potential and V the volume. In order to overcome this problem, in Ref. [80] a procedure 
was proposed in which the pressure is calculated from the derivative of the free energy per particle; 
in this way, the HVH theorem is automatically satisfied. The same drawback is present also in the 
relativistic DBHF. On the contrary, the HVH theorem is strictly fulfilled within the SCGF method 
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Figure 4: Free energy per particle for several values of the temperature T, for symmetric matter (upper 
panels) and pure neutron matter (lower panels). The V18 (left side) and SFHo (right side) EoSs are 
compared. See text for details. 


[198]. If only two-body forces are used, the results at the two-body correlation level in some cases 
are similar to the Brueckner ones, in some others they differ appreciably according to the forces used. 
The main difference with the Brueckner scheme is the introduction of the hole-hole propagation in the 
ladder summation, which gives a repulsive contribution. As a result, the critical temperature in the 
SCGF approach with Argonne Vis potential is found to be about T, ~ 11.6 MeV [198], depending on 
the adopted NN interaction. We remind that the values of the critical temperature depend both on the 
theoretical approach and on the NN interaction employed. For example, it turns out that the critical 
temperature within the DBHF scheme is definitely smaller than in the non-relativistic scheme, about 
10 MeV [199, 200]. This cannot be due to relativistic effects, since the critical density is about 1/3 of 
the saturation density, but to a different behaviour of the Dirac-Brueckner EoS at low density. This 
point remains to be clarified. However, there are experimental data from heavy-ion reactions that point 
towards a value of T. > 15 MeV [201]. For completeness, we also mention an empirical extension of a 
variational microscopic model based on the APR EoS [202]. 

The number of phenomenological EoSs available at finite temperature is also quite small. These are 
mostly based on density-functional theory (DFT), either on a relativistic version with various parame- 
terizations or a non-relativistic one based on Skyrme functionals. Those models are more extensively 
discussed in Sec. 3.5 below. 

In closing this section, we show in Fig. 4 a comparison between the microscopic BHF V18 EoS (left 
panels) and the general-purpose SFHo EoS (right panels). In particular we display the free energy per 
particle at different values of the temperature, for SNM (upper panels) and PNM (lower panels). In 
both cases we notice the typical Van der Waals behaviour of the isothermal curves, which are quite 
similar up to about twice the saturation density, but diverge at large density, being SFHo softer than 
V18. * Recently, EoSs at finite temperature incorporating hyperonic degrees of freedom have been also 


The original routine for the SFHo EOS is available at the CompOSE repository https: / /compose.obspm.fr 
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developed [78, 203-207], as well as those including a phase transition to quark matter [208-213]. 


3.5 General-purpose EoSs 


Besides the microscopic approaches, the so-called “general-purpose” EoSs are able to cover a wide range 
of densities, temperatures, and charge fractions, describing both clustered and homogeneous matter. 
These EoSs are therefore suitable for applications to SNe and mergers. However, at present, only a 
few of them are available and directly applicable to simulations; we list below the general-purpose EoSs 
with only nucleonic degrees of freedom currently used in astrophysical applications. 


e The Hillebrandt and Wolff (H&W) EoS [214] has been calculated using a Nuclear Statistical Equi- 
librium (NSE) network based on the model of [215] in the density range 10? — 3 x 102 g/cm?. At 
higher densities, the EoS is computed in the single-nucleus approximation, and the nuclear inter- 
action employed is the Skyrme interaction. This EoS is still used in recent numerical simulations 
[11]. 


e The Lattimer and Swesty (LS) EoS [216] is a very widely used EoS in numerical simulations. It 
models matter as a mixture of heavy nuclei treated in the single-nucleus approximation (SNA), 
a particles, free neutrons and protons, immersed in a uniform gas of leptons and photons. Nuclei 
are described within a medium-dependent liquid-drop model, and a simplified NN interaction of 
Skyrme type is employed for nucleons. With increasing density, shape deformations of nuclei 
(non-spherical nuclei and bubble phases) are taken into account by modifying the Coulomb and 
surface energies, and the transition to uniform matter is described by a Maxwell construction. 


e The Shen et al. (STOS) EoS [217] is also widely used. As the LS EoS, matter is described as a 
mixture of heavy nuclei (treated in SNA), a particles, and free neutrons and protons, immersed 
in a homogeneous lepton gas. For nucleons, a RMF model with the TM1 interaction [218] is 
used; o particles are described as an ideal Boltzmann gas with excluded-volume corrections. The 
properties of the heavy nucleus are determined by WS-cell calculations within the TF approach 
employing parameterized density distributions of nucleons and a particles. 


e The Furusawa EoS of [219] is based on a NSE model, including light and heavy nuclei. For nuclei, 
the liquid-drop model is employed, and the nuclear interaction is the RMF parameterization TM1 
[218]. This EoS has been applied in CCSN simulations to study the effect of light nuclei in 
[220], and the dependence of weak-interaction rates on the nuclear composition during stellar core 
collapse in [221]. 


e The SFHo and the SFHx EoSs [183] are based on the Hempel and Schaffner-Bielich (HS) [194] 
approach, which is inspired by the extended NSE model, that takes into account an ensemble of 
nuclei and interacting nucleons. Nuclei are described as classical Maxwell-Boltzmann particles, 
and nucleons are described within the RMF model employing new parameterizations fitted to some 
NS radius determinations. Binding energies are taken from experimental data or from theoretical 
nuclear mass tables. Excluded-volume effects are implemented in a thermodynamically consistent 
way so that it is possible to describe the transition to uniform matter. We remind that NSE models 
have been employed for typical conditions of core-collapse supernova [222-224]; for a comparison 
among the different methods, please refer to Ref. [225]. 


e The EoSs SHT [205] and SHO [226], are computed using different methods in different density- 
temperature domains. At high densities, uniform matter is described within a RMF model. For 
non-uniform matter at intermediate densities, calculations are performed in the (spherical) WS 
approximation, incorporating nuclear shell effects [227]. In this regime, matter is modelled as a 
mixture of one average nucleus and nucleons, but no a particles [227]. 
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The GRDF (generalized relativistic density functional) is based on a relativistic mean-field model 
of nuclear matter with density-dependent nucleon-meson couplings using the functional depen- 
dence introduced in Ref. [228]. Besides nucleons, electrons and muons with experimental masses, 
photons and nuclei are included as degrees of freedom. Two-nucleon correlations in the continuum 
are considered as effective resonances, and the dissolution of nuclei is described with the help of 
medium-dependent mass shifts. Masses of nuclei are taken from [229]. 


The EoS by Schneider, Roberts, and Ott (SRO) [230] is computed using the SLy4 parametrization 
[231]. The model includes nucleons, which are treated as non-relativistic particles; a particles, and 
photons, electrons and positrons, all treated as thermally equilibrated non-interacting relativis- 
tic gases. At low densities and temperatures nucleons may cluster into heavy nuclei computed 
within the SNA. Additional information on the open-source SRO EoS code can be found in 
https: / /stellarcollapse.org/SROEOS. 


The RG(SLy4) EoS corresponds to the extended NSE model proposed in Refs. [232, 233], where 
excluded volume effects between nuclear clusters and unbound nucleons are implemented via 
energy shifts of clusters binding energies. For nuclei for which experimental masses are known, 
the mass tables [229] are used. Then, up to the drip lines, evaluated masses model by Duflo and 
Zuker [234] are employed. Beyond drip lines, nuclear binding energies are described according to 
the liquid-drop-model-like parametrization of [235], corresponding to the SLy4 parameterization. 


The TNTYST EoS [195] is based on the variational many-body theory with realistic nuclear forces. 
For uniform matter, the EoS is constructed with the cluster variational method starting from the 
Argonne V18 two-body nuclear potential and the Urbana IX three-body nuclear potential. Non- 
uniform nuclear matter is treated in the Thomas-Fermi approximation. In the FT version, the 
variational approach is combined with the quantum approach for d, t, h and a, as well as the 
liquid-drop model for the other nuclei under the NSE assumption [219, 220]. 


The interested reader can refer to the CompOSE repository, https: / /compose.obspm.fr/, for addi- 


tional information regarding the practical use of the above listed EoS. Furthermore, we’d like to mention 
a few hybrid general purpose EoSs which are present in the repository, namely 
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e TheSTOS EoS has been used for describing the hadron phase, with the inclusion of a transition to 


quark matter [210] using a non-linear RMF model with the TM1 parametrization of the effective 
interaction and a bag model for the quark phase with three different values of the bag constant 
and the QCD coupling constant a,. The transition from the hadronic to the quark phase is 
done via a Gibbs construction. Non-uniform nuclear matter is calculated in the single-nucleus 
Thomas-Fermi approximation with parametrized density distributions in spherical Wigner-Seitz 
cells. Only neutrons, protons, alpha particles and a single heavy nucleus are considered. 


The relativistic density functional (RDF) formalism is used to describe homogeneous quark and 
hadron matter in a mean-field approximation [236], assuming a first-order phase transition from 
hadron to quark matter. The model was applied to NS configurations, and its current extension 
to finite temperature and arbitrary charge fractions has been employed to core-collapse supernova 
simulations [237] and binary NS merger simulations [238]. 


Structure and EoS of the neutron star crust 


Although the crust is only a small fraction of the star mass and radius, it plays an important role 
in various observed astrophysical phenomena such as pulsar glitches, quasiperiodic oscillations in soft 
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gamma-ray repeaters (SGR) and thermal relaxation in soft x-ray transients (SXT) [1, 239-243], which 
depend on the departure of the star from the picture of a homogeneous fluid. 

In modelling the NS crust, it is usually assumed that it has the structure of a regular lattice, in which 
nucleons can be either treated as a uniform system of interacting particles, or distributed within a defined 
shaped and sized cell. In the latter case, often the Wigner-Seitz (WS) approximation is used: matter 
is distributed in charge-neutral cells, which at lower densities are usually assumed spherical, centered 
around the positive-charged ion surrounded by a uniform electron and eventually free neutron gas, 
whereas at higher densities nuclei can be non-spherical and other geometries of the cell are considered. 
The minimization of the (free) energy of the system with respect to the variational variables, e.g. the 
nucleus atomic and mass number, the radius of the cell, and the free nucleon densities, under baryon 
number and charge conservation, allows the calculation of the EoS for the crust. If “pasta” phases 
are included, the minimisation is also performed on the shape of the cell. In the following we discuss 
the main features of the structure and crust EoS, and we refer the interested reader to Ref. [17] for a 
detailed treatment. 

Current ab-initio many-body calculations with realistic interactions are not affordable to describe 
clusterized matter, and therefore for the crust we have to rely on phenomenological models. As far 
as the outer crust is concerned, the classical way to determine the EoS is to use the so-called Baym- 
Pethick-Sutherland (BPS) model [3]. In this model, the outer crust is supposed to be made of fully 
ionised atoms arranged in a body-centered cubic lattice at T = 0 and to contain homogeneous crystalline 
structures made of one type of nuclides, coexisting with a degenerate electron gas. The EoS in each 
layer of pressure P is found by minimising the Gibbs free energy per nucleon, the only microscopic 
input being nuclear masses [1]. The BPS model has been updated by using modern nuclear data and 
theoretical mass tables [244], which were compared to check their differences concerning the neutron 
drip line, magic neutron numbers, the EoS and the sequence of neutron-rich nuclei up to the drip line 
in the outer crust. Pairing and deformation effects were also included. 

A very popular model to describe the inner crust was first introduced by Baym, Bethe and Pethick 
(BBP) [3], based on the compressible liquid drop model (CLDM) to take into account the effect of the 
dripped neutrons. Liquid-drop models were among the earliest to be used in astrophysical applications 
to treat non-uniform matter at zero and finite temperature, because of their applicability and limited 
computational efforts [3, 192, 245-251]). Those models parameterize the energy of the system in terms 
of global properties such as volume, asymmetry, surface, and Coulomb energy; their parameters are 
fitted phenomenologically. Moreover nucleons inside neutron-proton clusters and free neutrons outside 
are treated separately, and assumed to be uniformly distributed. Clusters have a sharp surface, and 
quantum-shell effects are neglected. A partially phenomenological approach, based on the CLDM, was 
developed by Lattimer and Swesty (LS) [216], who derived the EoS from a Skyrme nuclear effective 
force. Another EoS was developed by Shen et al. [217, 252] based on a nuclear RMF model. The crust 
was described in the Thomas-Fermi (TF) scheme using the variational method with trial profiles for 
the nucleon densities. The LS and Shen EoSs are widely used in astrophysical calculations for both NSs 
and supernova simulations due to their numerical simplicity and the large range of tabulated densities 
and temperatures. 

Douchin and Haensel (DH) [250] formulated a unified EoS for NS on the basis of the SLy4 Skyrme 
nuclear effective force [231], where some parameters of the Skyrme interaction were adjusted to re- 
produce the Wiringa et al. microscopic calculation of neutron matter [253] above saturation density. 
Hence, the DH EoS contains certain microscopic input. In the DH model the inner crust was treated in 
the CLDM approach. More recently, unified EoSs for NSs have been derived by the Brussels-Montreal 
group [254-257]. They are based on the BSk family of Skyrme nuclear effective forces [148]. Each force 
is fitted to the known masses of nuclei and adjusted among other constraints to reproduce a different 
microscopic EoS of neutron matter with different stiffness at high density. The inner crust is treated 
in the extended Thomas-Fermi approach with trial nucleon density profiles including perturbative shell 


26 


0.0005 - 
0.0004] —-— Shen-TM1 
- - - - BSk21 
"E —— BCPM 
E 0.0003 + 
u— 
> 
[0] 
=, 0.0002 - 
a 
0.0001 4 
0.0000 .0001 mn sg __—___- 
0.0000 0.0001 0.0002 0.0003 0.001 0.01 0.1 
-3 -3 
p, lim” p, lim” 


Figure 5: Pressure p as a function of the nucleon density for the outer crust (left) and the inner crust 
(right). Note that the x-axis scale on the right figure starts at p = 3 x 10—*fm"*, where the left figure 
ends. See text for details. Adapted from Ref. [258]. 


corrections for protons via the Strutinsky integral method. Analytical fits of these NS EoSs have been 
constructed in order to facilitate their inclusion in astrophysical simulations [257]. More recently, in 
[258] has been proposed a unified EoS (BCPM) from the outer crust to the core based on modern mi- 
croscopic calculations using the Argonne Vis potential plus TBFs computed with the Urbana model. To 
deal with the inhomogeneous structures of matter in the NS crust, they used a nuclear energy density 
functional that is directly based on the same microscopic calculations, and which is able to reproduce 
the ground-state properties of nuclei along the periodic table. The EoS of the outer crust requires the 
masses of neutron-rich nuclei, which are obtained through Hartree-Fock-Bogoliubov calculations with 
the new functional when they are unknown experimentally. To compute the inner crust, Thomas-Fermi 
calculations in Wigner-Seitz cells are performed with the same functional. We notice that the DH, 
BSk and BCPM approaches are the only ones which construct a unified theory able to describe on a 
microscopic level the complete structure of NS, from the outer crust to the inner core within the same 
theoretical approach. 

Quantal Hartree-Fock models have been also used for the NS crust, in which pairing and shell 
effects are properly taken into account [259-264]. Among those, we mention the Negele-Vautherin 
(NV) EoS as the mostly used one of this class. The main drawback of those calculations, which employ 
non-relativistic interactions, is that they are computationally very expensive. 

In Fig. 5, we display in the left panel an EoS sample for the outer crust, and in the right panel the 
one for the inner crust. The initial baryon density in the right panel corresponds to the last density 
shown in the left panel, where the EoS of the outer crust is plotted. The LS EoS [216], taken here in its 
Ska version (LS-Ska), and the Shen- TM1 EoS [217, 252] were computed with the Skyrme force Ska and 
the relativistic mean-field model TM1, respectively. In the two cases the calculations of masses are of 
semiclassical type and A and Z vary in a continuous way. Therefore, these EoSs do not present jumps 
at the densities associated with a change from an equilibrium nucleus to another in the composition. 
On the contrary those jumps are present for the BCPM and the BSk21 Skyrme nuclear effective force 
[148, 255-257] tabulated in [257]. The parameters of these forces were fitted to reproduce with high 
accuracy almost all known nuclear masses, and to various physical conditions including the neutron- 
matter EoS from microscopic calculations. We see that the BSk21 pressure is similar to the BCPM, with 
just some displacement around the densities where the composition changes from a nucleus to the next 
one. In the seminal work of BPS [3] the nuclear masses for the outer crust were provided by an early 
semi-empirical mass table. In Fig. 5 the corresponding EoS is seen to display a similar pattern with 
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the BCPM and BSk21 results. The LS-Ska EoS shows good agreement with the previously discussed 
models, with some departure from them in the transition to the inner crust. The largest discrepancies 
in Fig. 5 are observed with the Shen-TM1 EoS that in this region predicts a softer crustal pressure with 
density than the other models. 

The crustal pressure is an essential ingredient entering the Tolman-Oppenheimer-Volkoff (TOV) 
equations because of its implications for astrophysical phenomena such as pulsar glitches [243]. The 
pressure in the inner crust is provided by the free gases of the electrons and of the interacting dripped 
neutrons (aside from a correction from Coulomb exchange). We note, however, that the pressure 
obtained in a WS cell in the inner crust differs from the value in homogeneous npe matter in beta- 
equilibrium at the same average density owing to the lattice effects, which influence the electron and 
neutron gases. The lattice effects take into account the presence of nuclear structures in the crust and 
are automatically included in the self-consistent TF calculation. In the inner crust, the pressure from 
the BCPM functional is comparable in general to the results of the NV, BBP, DH-SLy4 and BSk21 
calculations. Particular agreement is observed in the inner crust regime between the BCPM and BSk21 
pressures up to relatively high crustal densities. In contrast, large differences are found when the BCPM 
pressure is compared with the values from LS-Ska and Shen-TM1 models. As the crust-core transition 
is approached, these differences can be as large as a factor of two, which may have an influence on the 
predictions of the mass-radius relationship of NSs, particularly in small-mass stars. 

All the approaches discussed above use the so-called single-nucleus approximation, i.e. the com- 
position of matter is assumed to be made of one representative heavy nucleus, usually the one most 
energetically favoured, along with alpha particles and unbound nucleons. It has been shown that this 
approximation has a small impact on thermodynamic quantities [265] with respect to the use of a full 
distribution of nuclei, though important differences might show up if the composition is dominated 
by light nuclei, or in the treatment of nuclear processes like electron captures in CCSNe. Indeed, the 
nucleus that is energetically favoured from thermodynamic arguments might not be the one with the 
highest reaction rate. 

At finite temperature, different configurations are expected and they can be calculated by using time- 
dependent Hartree-Fock models based on a wavelet representation [266, 267]. Further methods have 
also been developed, i.e. the classical and quantum molecular dynamics, where nucleons are represented 
respectively by point-like particles or wave packets [268-272], but these approaches are computationally 
very demanding. 

On the other hand, one can assume that the system is in thermodynamic equilibrium and use NSE 
models, where cluster degrees of freedom are introduced explicitly. These models suppose that the sys- 
tem is composed of a statistical ensemble of different nuclear species and nucleons in thermal, mechanical 
and chemical equilibrium. In the simplest version, NSE approaches treat the matter constituents as a 
mixture of non-interacting ideal gases, assuming they obey a Maxwell-Boltzmann statistics, although 
the Fermi-Dirac distribution is often adopted just for nucleons. NSE calculations require as input the 
nuclear binding energies, which can be either experimental, whenever available [273, 274], or theoret- 
ical (e.g. CLDM or EDF-based mass models). The main drawback of standard NSE-based models is 
that they neglect interactions and in-medium effects, that are essential for nuclear matter. Therefore, 
homogeneous matter cannot be correctly described by this kind of approaches. 

A possible way out is given by the extended NSE models, where the distribution of clusters is 
obtained self-consistently under conditions of statistical equilibrium, and interactions are properly in- 
cluded. For example, in-medium corrections of nuclear binding energies have been calculated for Skyrme 
interactions in [263, 275] within a local-density and extended TF approximation, respectively; also the 
screening of the Coulomb interaction due to the electron background is taken into account [223, 276]). 
The interactions between the cluster and the surrounding gas are often treated with an excluded-volume 
method (e.g. [220, 223, 277]). We notice that thermodynamic quantities are not very much affected 
by the presence of the ensemble of nuclei with respect to the single-nucleus approximation picture. 
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However, some quantitative differences arise in the matter composition, in particular concerning the 
contribution of light and intermediate mass nuclei (see, e.g. [232, 278]). Besides the first applications 
of a NSE model for the EoS of SN cores at low density, NSE models have been subsequently employed 
for conditions encountered in CCSN, e.g. in [194, 220, 222, 223, 232, 277, 279]. 


5 Structure and EoS of the core: Hyperons 


Since the pioneering work of Ambartsumyan and Saakyan [280], published in 1960, the presence of 
hyperons (baryons with a strangeness content) in the NS interior and their effects on the properties of 
these objects have been considered by many authors using either phenomenological (see e.g. Refs. [180, 
281-289]) or microscopic (see e.g. Refs. [290—302] approaches for the NS-matter EoS with hyperons. 
Contrary to terrestrial conditions, where hyperons are unstable and decay into nucleons through the 
weak interaction, the equilibrium conditions in NSs can make the inverse process happen. Hyperons 
may appear in the inner core of NSs at densities of about 2—3po. At such densities, the nucleon chemical 
potential is large enough to make the conversion of nucleons into hyperons energetically favourable. This 
conversion relieves the Fermi pressure exerted by the baryons and makes the EoS softer. Consequently, 
the mass of the star, and in particular its maximum value, is reduced. How much the EoS is softened, 
and how much the maximum mass is reduced depends on the attractive or repulsive character of the 
YN and YY interactions. In principle, attractive (repulsive) interactions cause an earlier (later) onset 
and larger (smaller) concentration of hyperons and, therefore, a stronger (more moderate) softening 
of the EoS and a larger (smaller) reduction of the maximum mass. However, it is well known (see 
e.g. Refs. [296, 297]) that, due to several compensation mechanisms, hyperons equalize the effect of 
different nucleonic interactions: a stiffer nucleonic EoS will lead to an earlier onset of hyperons, thus 
enhancing the softening due to their presence. Conversely, a later onset of a certain hyperon species 
will favour the appearance of other species leading also to a softer EoS. 

In the following we will review first some of the laboratory constraints of the nuclear EoS with 
hyperons derived from hypernuclear research. Then we will discuss the construction of YN and YY 
interactions on the basis of the meson-exchange and chiral effective field theories. Recent developments 
based on the so-called Vjows« approach and lattice quantum chromodynamics will also be addressed. 
Finally, we go over some of the effects of hyperons on the properties of neutron and PNSs with an 
emphasis on the so-called *hyperon puzzle". We discuss some of the solutions proposed to tackle this 
problem. We also re-examine the role of hyperons on the cooling properties of newly born NSs and on 
the development of the so-called r-mode instability. 


5.1 Laboratory constraints of the nuclear EoS with hyperons: hypernuclei 


The YN and YY interactions are still poorly constrained contrary to the NN one, which is fairly well 
known thanks to the large number of existing scattering data and known properties, such as masses and 
radii, of more than 3000 isotopes. Experimental difficulties associated with the short lifetime of hyperons 
in addition to the low-intensity beam fluxes have limited the number of AN and XN scattering events 
to just a few hundred [203—307] and that of EN events to very few. In the case of the YY interaction 
the situation is even worse because no scattering data exis,t at all. In view of this limited amount of 
data, which is not enough to fully constrain these interactions, alternative information on them can be 
obtained from the study of hypernuclei, bound systems composed of nucleons and one or more hyperons. 
In the following we present a brief review of the different production mechanisms of hypernuclei as well 
as a few aspects of their spectroscopy. 

Hypernuclei were first observed in 1952 [305], when a hyperfragment in a balloon-flown emulsion 
stack was found. Since then the use of high-energy accelerators as well as modern electronic counters 
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have lead to the identification of more than 40 single A-hypernuclei and a few double-A and single-= ones. 
The existence of single-X hypernuclei has not been experimentally confirmed yet without ambiguity, 
suggesting that the YN interaction is most likely repulsive. 

Single A-hypernuclei can be produced by several mechanisms such as (K~, m7) strangeness exchange 
reactions, where a neutron of the nucleus target hit by a K is changed into a A emitting a 7. The 
analysis of these reactions showed many of the hypernuclear characteristics such as, for example, the 
small spin-orbit strength of the YN interaction. The use of m? beams allowed (z+, K+) associated 
production reactions, where an ss pair is created from the vacuum and a K* and a A are produced 
in the final state. The electroproduction of hypernuclei by means of the reaction (e, e’A*) provides a 
high-precision tool for the study of hypernuclear spectroscopy due to the excellent spatial and energy 
resolution of the electron beams. Recently, the HypHI Collaboration at FAIR/GSI has proposed a new 
way to produce hypernuclei by using stable and unstable heavy ion beams [309]. The A and the ¿H and 
4H hypernuclei have been observed in a first experiment using a Li beam on a ^C target at 2 AGeV 
[310]. 

3. hypernuclei can also be produced in similar reactions. However, as we said, there is not yet an 
unambiguous experimental confirmation of their existence. The production of double-A hypernuclei, 
requires first to create a =~ through reactions like 


p+p> 8 +5? (36) 


or 
K na -K*. (37) 


Then it is necessary that the =~ be captured in an atomic orbit and interacts with the nuclear core, 
producing two A’s by hitting one of the protons of the nuclei 


E +poA+A. (38) 


This reaction releases about 30 MeV of energy that, in most of the cases, is equally shared between the 
two A’s, leading to the escape of one or both hyperons from the nucleus. Single-= hypernuclei can be 
produced by means of the reactions (36) and (37) and few of them have been in fact identified. The 
analysis of the experimental data from reactions such as *C(K”, K+)? Be [311] indicates an attractive 
=-nucleus interaction of about —14 MeV. We should mention here the recent observation [312] of a 
deeply bound state of the =~-'*N system with a binding energy of 3.87 + 0.21 MeV [313]. This event 
provides the first clear evidence of a deeply bound state of this system by an attractive =N interaction. 
Future = hypernuclei experiments are being planned at J-PARC. 

Double-strange hypernuclei are currently the best systems to investigate the properties of the baryon- 
baryon interaction in the strangeness S = —2 sector. The AA bound energy A Daa in double-A hyper- 
nuclei can be determined experimentally from the measurement of the binding energies of double and 
single-A hypernuclei as 

ABm = BMZ) —2Bi 2). (39) 


Emulsion experiments have reported the formation of a few double-A hypernuclei: $ He, A Be and 
1 B. From the subsequent analysis of these emulsion experiments, a quite large AA bound energy of 
around 4 to 5 MeV was deduced. We should also note that the identification of some of these double-A 
hypernuclei was ambiguous. Therefore, careful attention should be paid when using the data from this 
old analysis to put any kind of constraint on the AA interaction. However, a new H, He candidate with 
a AA bound energy ABra = 1.01 + 0.2101? MeV (recently corrected to ABaa = 0.67 + 0.17) was 
unambigously observed in 2001 at KEK [314]. 

Hypernuclei can be produced in excited states if a nucleon in a p or a higher shell is replaced 
by a hyperon. The energy of these excited states can be released either by emitting nucleons, or, 
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sometimes, when the hyperon moves to lower energy states, by the emission of y-rays. The detection 
of y-ray transitions in A hypernuclei has allowed the analysis of hypernuclear excited states with very 
good energy resolution. However, there have been some technical difficulties in the application of y-ray 
spectroscopy to hypernuclei mainly related with the detection efficiency of y-ray measurements and with 
the necessity of covering a large solid angle with y-ray detectors. The construction of large-acceptance 
germanium detectors, dedicated to hypernuclear y-ray spectroscopy, has allowed one to solve these 
issues somehow. There exist still, however, several weak points in hypernuclear y-ray spectroscopy. A 
number of single-particle A orbits are bound in heavy A hypernuclei with a potential depth of around 
30 MeV but the energy levels of many single-particle orbits are above the neutron and proton emission 
thresholds. Therefore, the observation of y-rays is limited to the low excitation region, maybe up to 
the A p-shell. The fact that the y-ray transition only measures the energy difference between two states 
is clearly another weak point, since single energy information is not enough to fully identify the two 
levels. The measurement of two y-rays in coincidence might help to resolve it. Systematic spectroscopic 
studies of single A hypernuclei indicate that the AN interaction is attractive [315]. 

Before closing this section, we would like to note that the YN and YY interactions can be tested 
experimentally by measuring the correlation in momentum space for Yp or YY pairs produced at 
colliders exploiting the so-called femtoscopy technique [316-323]. The correlation function between two 
hadrons is defined as the ratio of the distribution of relative momenta k* = |p; — p2|/2 for correlated 
and uncorrelated baryon pairs. In the case of an attractive interaction, the measured correlation will be 
larger than one, whereas if the interaction is repulsive or there exists a bound state, it will take values 
between zero and one. Theoretically, the correlation function can be expressed as [324, 325] 


C(k*) = / rs ua, (40) 


where S(r) is the distribution at the distance r at which particles are emitted (source) and W(k*, r) 
represents the relative wave function of the hadron pair of interest. The comparison of the measured 
correlation function with the theoretical one allows to test and improve the existing baryon-baryon 
potentials. 


5.2 The hyperon puzzle 


Although, energetically, the presence of hyperons in NSs seems to be unavoidable, the softening that 
their presence induces on the EoS leads to maximum masses which are incompatible with observa- 
tions. The solution of this problem, commonly known in the literature as the “hyperon puzzle,” is 
not easy. It has become a subject of very active research (see e.g. Ref. [326] and references therein for 
a comprehensive review), since the measurements of unusually high masses of the millisecond pulsars 
PSR J1903+0327 (1.667 3: 0.021 Mọ) [327], PSR J1614-2230 (1.928+0.017 Mọ) [328], PSR J0348+0432 
(2.01 + 0.04 Mo) [329] and the most recent PSR J0740--6620 (2.14*019 Mo) [330], which rule out al- 
most all currently proposed EoSs with hyperons (both microscopic and phenomenological). To solve 
this problem, a mechanism is necessary that could provide the additional repulsion needed to make the 
EoS stiffer and therefore the maximum mass compatible with the current observational limits. Three 
of the mechanisms proposed that could provide such additional repulsion are: (i) the inclusion of a re- 
pulsive YY interaction through the exchange of vector mesons [331-334], (ii) the inclusion of repulsive 
hyperonic TBFs [335-342], or (iii) the possibility of a chiral phase transition and/or a phase transition 
within color superconducting states at densities below the hyperon threshold [343-350]. The possible 
appearance of other hadronic degrees of freedom, such as the A isobar or meson condensates, that can 
push the onset of hyperons to higher densities has also been considered. An interesting possibility to 
circumvent the problem is the so-called two-families scenario proposed in [351, 352], in which stars made 
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of hadrons are stable only up to (1.5— 1.6) Mo, while most massive compact stars are entirely made of 
strange quark matter. In the following, we briefly review some of these possible solutions. 


5.2.1 Hyperon-hyperon repulsion 


This solution, mainly explored in the context of RMF models (see e.g. Refs. [331-334]), is based on 
the known fact that in meson-exchange models of nuclear forces, vector mesons are responsible for 
the repulsion felt by nucleons at short distances, while the o meson generates the intermediate-range 
attraction. Therefore, it has been argued that if the interaction of two hyperons, mediated through 
the exchange of a vector meson, is repulsive enough or the attraction driven by the e meson is weak 
enough, then it can be possible to reconcile the current high pulsar mass observations with the existence 
of hyperons in NSs. The existence of hypernuclei, however, shows that at least the AN interaction is 
attractive [315]. Consequently, to be consistent with hypernuclear data, repulsion in these models is 
included only through the exchange of the hidden strangeness @ vector meson that couples only to 
the hyperons. In such a way, the onset of hyperons is shifted to higher densities, and it is possible to 
obtain NSs with maximum masses larger than 2 M; and with a non-negligible fraction of hyperons in 
their interior. A question, however, immediately arises: how does this additional YY repulsion affect 
the binding energy of two As in double-A hypernuclei ? Although the answer is not simple, it has 
been shown recently [353] that with currently available hypernuclear experimental data and due to the 
lack of stronger constraints on the asymmetric nuclear-matter EoS, in RMF models it is still possible 
to find a wide range of values of the Ad coupling that predict a AA binding energy compatible with 
that derived from experimental data on H, He, and that simultaneously lead to the prediction of NS 
maximum masses compatible with the 2 Mo constraint. The interested reader is referred to any of the 
works (see e.g. Ref. [326] and references therein) that have explored this solution over the last years for 
further information. 


5.2.2 Hyperonic three-body forces 


It is well known from conventional nuclear physics that TBFs are fundamental to reproduce properly 
the properties of few-nucleon systems as well as the empirical saturation point of SNM. Then, it is 
natural to think that TBFs involving one or more hyperons (i.e. NNY, NYY and YYY) could also play 
an important role in the determination of the NS-matter EoS and contribute to the solution of the 
hyperon puzzle. Indeed if hyperonic TBFs are repulsive enough then they can make the EoS stiffer 
at higher densities and, therefore, make the maximum mass of the star compatible with observations. 
This solution was suggested years ago before the observation of NSs with ~ 2 Mọ (335, 336] and it has 
been considered later by other authors [337-342]. A general consensus regarding the role of these forces 
on the solution of the hyperon puzzle, however, does not exist yet. In Refs. [338, 339], for instance, 
it was found that the inclusion of hyperonic TBFs is enough to obtain hyperon stars with ~ 2 Mo. 
However, in Ref. [337] it was shown that the largest maximum mass that they can support is ~ 1.6 Mo. 
Finally, the results of Ref. [340] were not conclusive enough because they depend strongly on the ANN 
interaction employed. It seems, therefore, that hyperonic TBFs are not the full solution of the hyperon 
puzzle, but they can probably contribute to it in an important way. We want to mention here that the 
effect of NNY forces, derived in the context of chiral effective field theory, on the properties of hyperonic 
matter has been recently studied in [341] and [342]. The interested reader is referred to all the works 
quoted, Refs. [335-342], for the specific details of the calculations. 


5.2.3 Quarks in neutron stars 


It has been suggested that an early chiral phase transition and/or a phase transition within color 
superconducting states at densities below the hyperon threshold could provide a solution to the hyperon 
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puzzle. Compact stars which possess a quark-matter core, either as a mixed phase of deconfined quarks 
and hadrons, or as a pure quark matter phase, surrounded by a shell of hadronic matter are called 
hybrid stars. Massive stars could therefore actually be hybrid stars with a stiff quark-matter core. The 
question that arises then is whether quarks can provide sufficient repulsion to support a 2 M; NS. A 
brief summary of some important conclusions drawn in the last years using the massive NS constraint 
is given in the following. This summary does not pretend to be exhaustive and therefore we apologize 
in advance to those groups whose results are not mentioned here. 

Following the discovery of the massive NS EXO 0748-676, it was argued [343] that hybrid or strange 
stars (i.e. compact stars made of deconfined u, d, s quark matter) can reach a mass of 2 Mo, as formerly 
demonstrated within the framework of the MIT bag model, perturbative QCD as well as the Nambu- 
Jona-Lasinio (NJL) model. A systematic study of the consequences of the NS mass limit on the 
properties of quark stars and hybrid stars with a pure quark core was performed in [344]. Using 
an extended quark bag model, the authors concluded that massive strange stars require strong QCD 
corrections and large contribution from color superconductivity. For the case of hybrid stars, the EoSs 
with quark-hadron phase transition are compatible with the 2 Mo mass constraint, provided the quarks 
are strongly interacting and in a color-superconducting phase [344, 345]. Ref. [346] defined a hadronic 
flavor-SU(3) model to combine hadronic and quark degrees of freedom in a unified way. They obtained 
a maximum mass of 2.06 Mọ for models with hyperons only, and reported a mass reduction of 10% 
on inclusion of quark degrees of freedom. A further increase of about 20% in the maximum mass 
was achieved by considering rotation. Ref. [347] used the heavy NS mass observation to put general 
constraints on the two-flavor color-superconducting (2SC) and color-flavor-locked (CFL) phases of quark 
matter in the core of hybrid stars. They deduced that for thermodynamic stability, both a stiff hyperon 
repulsion at high baryon density accompanied by a stiff quark matter EoS are required to reproduce this 
result. Ref. [348] described a large number of possible parametrizations based on the NJL model that 
can reach masses up to 2 M¿. They demonstrated that to reconcile NS observations and heavy-ion flow 
data, large values of the diquark and scalar couplings are necessary, as well as repulsion arising from 
the vector-meson interaction. Ref. [349] succeeded in constructing stable NS configurations with mass 
> 2 Mo using an EoS exploiting features of phenomenological RMF density functionals at low densities, 
and the NJL model of superconducting quark matter at high densities. The constraints were satisfied 
subject to the condition that the nuclear EoS at post-saturation density is sufficiently stiff, followed 
by a transition to quark matter at a few times saturation density, with substantial repulsive vector 
interactions in quark matter. Ref. [350] studied the possible appearance of hyperons and strange quark 
matter in NSs subject to the constraints of observations of heavy compact stars and flow constraints from 
heavy ion collisions (HICs). Applying a color superconducting three-flavor NJL model for the quark 
sector, and a density-dependent hadron field model and the DBHF theory for the hadronic sector, they 
showed that it is possible to have deconfined quark matter in the core of massive stars. 


5.2.4 A isobar and meson condensation in neutron stars 


The appearance of other hadronic degrees of freedom such as for instance the A isobar or meson 
condensates that can push the onset of hyperons to higher densities has also been considered as a 
possible solution to the hyperon puzzle. The presence of the A isobar in NSs has been usually neglected 
because its threshold density was found to be larger than the typical densities prevalent in the NS core. 
However, it has been recently shown [354] that the onset density of the A isobar depends crucially on 
the slope parameter L of the nuclear symmetry energy. Using recent experimental constraints and a 
state-of-the-art EoS, these authors have shown that the A isobar could actually appear at densities 
below the hyperon threshold. However, they found that as soon as the A appears, the EoS becomes 
also considerably soft and consequently the maximum mass is reduced to values below the current 
observational limit, giving rise to what they have named “the A puzzle.” The possibility that pions or 


33 


—==="=—3 
T=0 MeV S/A =1 S/A = 2 


0.14 Be este ET | fe oe eee 


0.14 af ai a d ge A eet Us 
d Å e SE GE, LN 


^ iy KEN r 
n. Fi PE Ei LS ` 
D 8 M y sg å Kei 
D 1 = $ T. E 
0.014 j s ; PU E 
` Te. ; ` 
H SJ ! 
b d "Al 
0.001 Lt —— ty : å : r 
00 02 04 06 08 00 02 04 06 08 00 02 04 06 08 1.0 


p [fm?] 


T T T T T T T T T T T T 
dl El m | 


T=0 MeV S/A=1 S/A = 2 
0.1 LE MM CL A E 
DLL DID f cL E 
ss ZZ L Ó 1 1 d 
e " "P d 
0.01 4 ---- v, PE Pa 
=.-. H š = m 
eX 1 y E 
O TT A å s ` i 
x^ 0.001 r : —— — —— r — —— T 


0.14 =. 


0.01 : : 
: d pro 
i ig For 
ni PN Rå B le 
0.001 r : —— r mr r —— r 
0.0 02 04 06 08 00 02 04 06 08 00 02 04 06 08 1.0 
3 
p [fm] 


Figure 6: Composition for neutrino-free (upper plot) and neutrino-trapped (lower plot) matter for 
several fixed values of the entropy per particle S/A — 0,1,2 with (lower rows) and without (upper 
rows) hyperons. Figure adapted from Ref. [361]. 


kaons form Bose-Einstein condensates in the inner core of NSs has also been extensively considered in 
the literature [355-360]. However, as in the case of the hyperons or the A isobar, the appearance of a 
meson condensate induces also a strong softening of the EoS and therefore leads to a reduction of the 
maximum NS mass to values also below the current observational limits. 


5.3 Effect of hyperons on proto-NSs and neutron star cooling 


Thermal effects and neutrino trapping characterise a newly born proto-NS (PNS). If hyperons are 
present, their number is on average smaller and their onset is shifted to higher densities than in cold 
NSs. Consequently, the EoS is stiffer in comparison with the cool and neutrino-free case. All these 
features are illustrated in Figs. 6 and 7 (adapted from Ref. [361]), where the composition and Foi of 
neutrino-free and neutrino-trapped (Y, — 0.4) matter with and without hyperons are shown for several 
fixed values of the specific entropy. 
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Figure 7: EoS for neutrino-free and neutrino-trapped matter for several fixed values of the entropy per 
particle S/A — 0,1,2 with and without hyperons. Figure adapted from Ref. [361]. 


An important consequence of neutrino trapping in NSs with hyperons is the possible delayed collapse 
of the hot newly born NS to a black hole, as was already pointed out in [363]. This possibility is 
illustrated in Fig. 8, which shows the gravitational mass Mg of the star as a function of its corresponding 
baryonic mass Mg for the BHF calculation of Ref. [562]. When hyperons are present in the star (left 
panel), the deleptonization lowers the range of gravitational masses that can be supported by the EoS 
from 1.59 Mo to 1.28 Mo (see dotted horizontal lines in the figure). The NS baryonic mass can be 
considered to a good approximation constant during the evolution from the initial PNS configuration 
to the final neutrino-free one, because most of the matter accretion on the forming NS happens in the 
very early stages after its birth (t S 1s). For the particular model calculation shown in the figure, 
PNSs which at birth have a gravitational mass between 1.28 Mo and 1.59 Mo (a baryonic mass between 
1.40 Mo and 1.72 Mọ) will be stabilized by neutrino effects long enough to carry out the nucleosynthesis 
accompanying a type-II supernova explosion. Once neutrinos have left the star, the EoS softens and it 
cannot any more support the star against its own gravity. Thus the newborn NS collapses into a black 
hole [12, 364, 365]. Conversely, when nucleons are considered to be the only relevant baryonic degrees of 
freedom present in the NS core, a black hole is unlikely to be formed during the deleptonization since in 
this case (see right panel) the gravitational mass increases during this stage which happens at (almost) 
constant baryonic mass. If a black hole were to form from a star with only nucleons, it is much more 
likely to form during the post-bounce accretion stage. 

The presence of hyperons may also affect the cooling of NSs, since it provides fast additional cooling 
mechanisms, such as for example direct and modified hyperonic Urca processes, 


Y>B+l4+y, B+l5Y+y, (41) 
B'+Y+3B+B+1+7, B+B+l>B'+Y+v. (42) 


Such additional cooling mechanisms, however, can lead to predicted surface temperatures much lower 
than those observed (see Sec. 7), unless they are suppressed by hyperon pairing gaps. Therefore, 
the study of hyperon superfluidity becomes of particular interest since it could play a key role in the 
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Figure 8: Gravitational mass as a function of the baryonic mass for neutrino-free (solid lines) and 
neutrino-trapped (dashed lines) matter. Panel (a) shows the results for matter including nucleons and 
hyperons, whereas panel (b) shows the results for pure nucleonic matter. Dotted horizontal and vertical 
lines show the window of metastability in the gravitational and baryonic masses. Figures adapted from 
Ref. [362]. 


thermal history of NSs. However, whereas the presence of superfluid neutrons in the inner crust of 
NSs and superfluid neutrons together with superconducting protons in their quantum fluid interior is 
well established and has been the subject of many studies (Sec. 7.7), a quantitative estimation of the 
hyperon pairing gaps has not received much attention and just a few calculations exist in the literature 
[366-372]. 


5.4 Hyperons and the r-mode instability of NSs 


It is well known that the Kepler frequency Qx is the absolute maximum rotational frequency of a NS, 
above which matter is ejected from the star’s equator [373, 374]. Instabilities caused by several types of 
perturbations may prevent NSs from reaching rotational frequencies as high as Qx, imposing therefore 
more stringent limits on their rotation [375]. Among the different instabilities that can operate in a 
NS, the so-called r-mode instability [376, 377], a toroidal mode of oscillation whose restoring force is 
the Coriolis force, is particularly interesting. The r-mode instability leads to the emission of GWs 
in hot and rapidly rotating NSs through the Chandrasekhar-Friedman-Schutz mechanism [378-381]. 
Gravitational radiation makes an r-mode grow, while viscosity stabilizes it. Therefore, if the driving 
time of the gravitational radiation is shorter than the damping time associated to viscous processes, 
the r-mode becomes unstable. In this case, a rapidly rotating NS could transfer a significant fraction 
of its angular momentum and rotational energy to the emitted GWs, whose detection could provide 
invaluable information on the internal structure of the star and constraints on the EoS. 

The main dissipation mechanisms of r- and other pulsation modes of NSs are the bulk (£) and shear 
(n) viscosities. Bulk viscosity is the dominant one at high temperatures (T > 10%K) and therefore 
it is important for hot young NSs. It is the result of the variations in pressure and density induced 
by the pulsation mode, variations that drive the star away from beta-equilibrium. Energy is then 
dissipated as the weak interaction tries to reestablish the equilibrium. In the absence of hyperons or 
other exotic components, the bulk viscosity of NS matter is mainly determined by the reactions of 
direct and modified Urca processes. However, as soon as hyperons appear, new mechanisms such as 
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Figure 9: Panel (a): r-mode instability region for a pure nucleonic and a hyperonic star with 1.27 Mo. 
The frequency of the mode is taken as w = 10*s^!. Panel (b): Bulk viscosity as a function of the 


density for T = 10? K and w = 10*s^!. Contributions of direct and modified nucleonic Urca processes 
as well as from the weak non-leptonic process n +n + p+ M^ are included. 


weak non-leptonic hyperon reactions 
N+N>-N+Y,N+YeY+Y, (43) 
strong interactions 
Yo Yaa Yo N+438Y+Y, Y+Y8ƏY+Y, (44) 


or direct and modified hyperonic Urca reactions, Eqs. (41) and (42), contribute to the bulk viscosity and 
dominate it for densities p > 2 — 3p9. Hyperon bulk viscosity has been considered by several authors, 
see e.g. Refs. [382-395]. 

The time dependence of an r-mode oscillation is given by e ) being w the frequency of the 
mode and 7(Q, T) an overall time scale describing both the exponential growth of the mode due to GW 
emission as well as its decay due to viscous damping. It can be written as 


iwt—t/7(Q,T 


1 1 1 1 


HD RO a nor’ i9) 


where Taw, Te and Ty are the time scales associated to the GW emission and the bulk and shear viscosity 
dampings, respectively. If raw is shorter than both ze and 7,, the mode will exponentially grow, while 
in the opposite case it will be quickly damped away. For each star at a given temperature T' one can 
define a critical angular velocity Q. as the smallest root of the equation 
1 
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This equation defines the boundary of the so-called r-mode instability region. A star will be stable 
against the r-mode instability if its angular velocity is smaller than its corresponding Qe. On the 
contrary, a star with Q > Q, will develop an instability that causes a rapid loss of angular momentum 
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through gravitational radiation until its angular velocity falls below the critical value. As an example, 
in panel (a) of Fig. 9 is presented the r-mode instability region for a pure nucleonic (black solid line) and 
a hyperonic (red dashed line) star with 1.27 M¿. The contributions to the bulk viscosity from direct 
and modified nucleonic Urca processes as well as from the weak non-leptonic process n +n €» p+ M^ 
included in the calculation are shown in panel (b) of the figure. Clearly the r-mode instability is smaller 
for the hyperonic star, because of the increase of the bulk viscosity due to the presence of hyperons, 
which makes the damping of the mode more efficient. 


6 Constraints on the EoS 


6.1 Nuclear experiment constraints of the EoS 


It is well known that around saturation density po and isospin asymmetry 8 = (~n—Pp)/(Pn+Pp) = 0, the 
nuclear EoS can be characterized by a set of a few isoscalar (Eo, Ko, Qo) and isovector (S, L, Ksym; Qua) 
parameters. Ep is the energy per particle of SNM at po, Ko the incompressibility parameter, Qo the 
so-called skewness, So the value of the nuclear symmetry energy at po, L the slope of the symmetry 
energy, Ksym the symmetry incompressibility and Qsym the third derivative of the symmetry energy with 
respect to the density. These parameters can be constrained by nuclear experiments and are related to 
the coefficients of a Taylor expansion of the energy per particle of asymmetric nuclear matter in density 
and isospin asymmetry, 


E 1 1 1 1 
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where x = (p — po)/3po- 


6.1.1 Constraints of the isoscalar parameters 


Measurements of density distributions [396] and nuclear masses [273] yield po = 0.15 — 0.16 fm? and 
Ey = —16 + 1 MeV, respectively. Although its extraction is complicated and not unambiguous, the 
value of Ko can be determined experimentally from the analysis of isoscalar giant monopole resonances 
in heavy nuclei. Results of Ref. [397], for instance, suggest Ko = 240 + 10 MeV, whereas in Ref. [398] a 
value of K = 24848 MeV was reported, and Ref. [399] gives K = 210430 MeV. Recently, in Ref. [400] it 
was shown that the third derivative M of the energy density of SNM is constrained by giant monopole 
resonance measurements not at saturation density but rather around what has been called crossing 
density p = 0.11 fm ?. The authors of this work found M = 1100 + 70 MeV, which, once extrapolated 
to po, implies Kg = 230 + 40MeV. We note that the value of the skewness parameter Qo is more 
uncertain and not very well constrained yet, being its estimate in the range —500 < Qo < 300 MeV. 

A rather “soft” EoS, i.e. a lower value of Ko [401], would be pointed out by HIC experiments. 
Constraints on the nuclear EoS for SNM from HIC data of the KaoS experiment [402] and flow data 
[403] are shown in the left panels of Fig. 10, together with the predictions of some microscopic (solid 
lines in the upper panels) and phenomenological (broken lines in the lower panels) EoSs considered. 
Any reliable model for the EoS should pass through the region defined by the experimental data. We 
notice that most of the EoSs considered are in general compatible with these data, with some exception 
in the case of the phenomenological ones, which are too repulsive and some of the microscopic ones 
(BOB, DBHF, V18), which are only marginally compatible. We should point our that the constraints 
inferred from HICs are however model dependent, since the analysis of the measured data requires the 
use of transport models. For completeness, in the right panels of the figure we show the astrophysical 
constraints (see next section) on the EoS for beta-stable matter imposed by the GW170817 event [8]. 
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Figure 10: The pressure as a function of the baryon density for symmetric (left panels) and beta-stable 
matter (right panels). Results for microscopic (phenomenological) models are shown in the upper 
(lower) panels. Constraints from HIC data of the KaoS experiment (orange bands) and flow data (gray 
bands) are shown together with the limits deduced by the GW170817 event (blue bands in the right 
panels). 


6.1.2 Constraints of the isovector parameters 


The isovector parameters of the nuclear EoS can be constrained experimentally from isospin diffusion 
measurements [405], the analysis of giant [406] and pygmy [407, 408] resonances, isoscaling [409], isobaric 
analog states [235], pion [410] and kaon [411] production in HICs or measurements of the neutron skin 
thickness in heavy nuclei [162, 412-418]. Astrophysical observations can also be used to constrain these 
parameters. It has been shown, for instance, that the slope parameter L of the symmetry energy is 
correlated with the radius [419] and the tidal deformability [420, 421] of a 1.4 Mg NS, and that precise 
and independent measurements of the radius and the tidal deformability from multiple observables of 
NSs can potentially pin down the correlation between Ksym and L and thus the high-density behaviour of 
the nuclear symmetry energy [122]. Nevertheless, whereas Sy is more or less well established ( 30 MeV), 
the values of L, and specially those of Ksym and Qsym, are still uncertain and poorly constrained. 
Combining different data, for instance, the authors of Ref. [423] give 29.0 < So < 32.7MeV and 
40.5 < L < 61.9 MeV, while in Ref. [187] it was suggested 30.2 < So < 33.7 MeV and 35 < L < 70 MeV. 
Why the isovector part of the nuclear EoS is so uncertain is still an open question whose answer is related 
to our limited knowledge of the nuclear force and in particular to its spin and isospin dependence. 

We show in Fig. 11 the slope of the symmetry energy L as a function of the symmetry energy at 
saturation So. The different symbols show the predictions from microscopic approaches (black circles), 
the Skyrme EoS (green triangles) and the NLWM (red squares) and DDM models (blue diamonds). 
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Figure 11: Slope of the symmetry energy L as a function of the symmetry energy at saturation So. 
Full symbols show the predictions from microscopic approaches (black circles), the Skyrme EoS (green 
triangles) and the NLWM (red squares) and DDM models (blue diamonds). The shaded areas represent 
experimental bands whereas the dashed line shows the unitary limit constraint determined in Ref. [404]. 


Currently available experimental data from neutron-skin thickness in Sn isotopes [124], isospin diffusion 
in HICs [125], electric dipole polarizability [126], isobaric-analog-state (IAS) phenomenology combined 
with the skin-width data [187] or the Bayesian analysis of mass and radius measurements in NSs [183] 
are also indicated in the figure. The dashed line shows the unitary limit constraint determined in 
Ref. [104], which combined with the experimental and observational data indicates that only values of 
(So, L) to the right of this line are permitted. As can be seen, all the microscopic and most of the 
phenomenological models considered fulfill all these constraints. 


6.2 Astrophysical constraints 


The observation of NSs provides the main astrophysical constraints on the nuclear EoS. An enormous 
amount of data on different NS observables has been collected after fifty years of observations. From 
these observables it is possible to infer valuable information on the internal structure of these objects and 
therefore also on the nuclear EoS, which is the only ingredient needed to solve the structure equations 
of NSs. In the following we shortly review some of these observables. 


6.2.1 Masses 


NS masses can be inferred directly from observations of binary systems and likely also from supernova 
explosions. In any binary system, there exist five orbital parameters (usually known as Keplerian 
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parameters), which can be precisely measured. These are: the eccentricity of the orbit e, the orbital 
period P,, the projection of the pulsar’s semi-major axis on the line of sight x = a,sini/c, where i 
is the inclination of the orbit, and the time Tọ and longitude wo of the periastron. With the help of 
Kepler’s third law, these parameters can be related to the masses of the NS (M,) and its companion 
(M.) through the so-called mass function 


(Mesini ^ Bag 
MIMI? 


TOM, Me, i) = (48) 


where v, = 2raısini/P, is the projection of the orbital velocity of the NS along the line of sight. 
One cannot proceed further than Eq. (48) if only one mass function can be measured for a binary 
system, unless additional assumptions are made. Deviations from the Keplerian orbit due to general 
relativity effects can be detected fortunately. These relativistic corrections are parametrized in terms 
of one or more post-Keplerian parameters. The most significant ones are: the advance of the periastron 
of the orbit (w), the combined effect of variations in the transverse Doppler shift and gravitational 
redshift around an elliptical orbit (Y), the orbital decay due to the emission of quadrupole gravitational 
radiation (P,), and the range (r) and shape (s) parameters that characterize the Shapiro time delay of 
the pulsar signal as it propagates through the gravitational field of its companion. These post-Keplerian 
parameters can be written in terms of measured quantities and the masses of the star and its companion 
as [427]: 
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where n = 27/P, is the orbital angular frequency and Ta = GM, = 4.925490947 x 10 9s. The 
measurement of any two of these post-Keplerian parameters together with the mass function f is 
sufficient to determine uniquely the masses of the two components of the system. An example of a 
high-precision mass measurement is that of the famous Hulse-Taylor binary pulsar [428] with measured 
masses M, = 1.4408 + 0.0003 Mo and M. = 1.3873 + 0.0003 Mo. Other examples are those of the 
millisecond pulsars PSR. J1614-2230 [328], PSR. J0348+0432 [329] and the most recently observed PSR 
J0740+6620 [330] with masses M, = 1.928 + 0.017 Mo, M, = 2.01 + 0.04 Mo and M, = 2.14*049 Mo, 
respectively. These are binary systems formed by a NS and a white dwarf. The measurement of these 
unusually high NS masses constitutes nowadays one of the most stringent astrophysical constraints on 
the nuclear EoS. 


6.2.2 Radii 


NS radii are very difficult to measure, mainly because NSs are very small objects and are very far 
away from us (e.g. the closest NS is probably the object RX J1856.5-3754 in the constellation Corona 
Australis, which is about 400 light-years from earth). Direct measurements of radii do not exist. 
However, a possible way to determine them is to use the thermal emission of low-mass x-ray binaries. 
The observed x-ray flux F and temperature T, assumed to be originated from a uniform black body, 
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together with a determination of the distance D of the star can be used to obtained the radius through 


the following implicit relation 
F D? 2GM 
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where c is the Stefan-Boltzmann constant and M the mass of the star. The major uncertainties in the 
measurement of the radius through Eq. (54) originate from the determination of the temperature, which 
requires the assumption of an atmospheric model, and the estimation of the distance of the star. The 
analysis of present observations from quiescent low-mass x-ray binaries is still controversial. Whereas 
the analysis of Ref. [183, 429] indicates NS radii in the range of 10.4 — 12.9 km, that of Ref. [430, 431] 
points towards smaller radii of ~ 10 km or less. 

We would like to note here that the simultaneous measurement of both mass and radius of the 
same NS would provide the most definite observational constraint on the nuclear EoS. Very recently the 
NICER mission has reported a Bayesian parameter estimation of the mass and equatorial radius of the 
millisecond pulsar PSR J0030+0451 [432, 433]. The values inferred from the analysis of the collected 
data are, respectively, 1.347912 Mo and 12.71*11$ km. We want to stress here that the measurement 
performed by the NICER mission does not make use of Eq. (54) and constitutes the first model- 
independent one, since only the geometry of the hot spots of the NS and general relativity effects enter 
in the determination of the star radius. 


6.2.3 Gravitational waves 


On August 17 2017, a GW signal from a binary NS merger was detected for the first time by the 
Advanced LIGO and Advanced VIRGO collaborations [8], inaugurating with this event (known as 
GW170817) a new era in the observation of NSs. GWs originated during the coalescence of two NSs 
or a black hole and a NS constitute now a new and valuable source of information on the EoS and 
internal structure of NSs. In particular, the so-called tidal deformability A, or equivalently the tidal 
Love number ka of a NS [134—436], can provide priceless information and constraints on the related EoS, 
because it depends strongly on the compactness = M/R of the object. The Love number, defined as 
à A = 3 GA 88% 
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with A = A/M? and 
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F = 68(2 — yn) + 68*(5yg — 8) + 48 (13 — 11yr) + 48 (3yr — 2) + 88 (1 + yn) + 3zIn(1— 28), 
can be obtained by solving the TOV equations 
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together with the additional equation [437] 
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We note that yr = y(R) is appearing in the definitions of z and F given above. Equations (56), (57) 
and (58) are solved using the EoS p(e) as the only input needed, and imposing the following set of 
boundary conditions 


[p, m, y] (r = 0) = [Pes 0, 2] : (60) 


The average tidal deformability of an asymmetric binary NS system with mass asymmetry q = 
M2/M and known chirp mass M, + (M,M3)P/(M, + M3)!/^, characterizing the GW signal waveform, 
is given by 
16 (1 + 129)A1 + (q + 12)A2 
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A value of M. = 1.186 i Bn: Mo was obtained from the analysis of the GW170817 event [8, 438, 439], 


corresponding to Mı = Ma = 1.365 M; for a symmetric binary system, q = 0.73 — 1 and A « 730 from 
the phase-shift analysis of the observed signal. We stress that the GW170817 observation put a strong 
constraint on the radius of a 1.4 Mo NS. Requiring both NSs to have the same EoS, this leads to the 
constraints 70 < Ay4 < 580 and 10.5 < R4 < 13.3km [438] for this radius. In Ref. [440] a general 
polytropic parametrization of the EoS compatible with perturbative QCD at very high density was 
used, and the constraint A44 < 800 yielded a similar upper limit Rj 4 < 13.4 km. 

The high luminosity of the kilonova AT2017gfo following the NS merger event imposes also a lower 
limit on the average tidal deformability, A > 400, which was deduced in order to justify the amount 
of ejected material being heavier than 0.05 M¿. This lower limit, which was used in Refs. [141-444] in 
order to constrain the EoS, could indicate that Rı4 Z 12km. However, this constraint of Á has to be 
taken with great care and, in fact, it has been recently revised to A > 300 [445], but considered of limited 
significance in Ref. [146]. We thus notice that the determination of the average tidal deformability of 
the binary NS system GW170817 has imposed constraints on the NS radii, to lie between about 12 
and 13 kilometers |447]. 'This is complementary to the mass-radius measurement by NICER mentioned 
before, and contributes to selecting the suitable EoS. 

Before discussing other astrophysical constraints, we show in Fig. 12 the mass-radius relations ob- 
tained with different microscopic (solid lines) and phenomenological (broken lines) EoSs. Most of the 
models considered, except the soft microscopic UIX and FSS2GC, predict values of the maximum mass 
larger than 2 Mo, therefore being compatible with current observational limits [328—330]. The mass of 
the most heavy pulsar PSR J07404-6620 [330] observed until now is also shown in the figure, together 
with the constraints of the NICER mission [132, 433] mentioned above and those from the GW170817 
event [8]. Some recent theoretical analyses of this event indicate an upper limit of the maximum mass of 
~ 2.2 — 2.3 Mo [148-450], with which several of the microscopic and phenomenological EoSs considered 
would be compatible as well. 


6.2.4 Rotational periods 


The majority of the known NSs are observed as radio pulsars which exhibit in most of the cases very 
stable rotational periods. Highly accurate pulsar timing has allowed observers to measure the rotational 
period P and the first time derivative P and sometimes even the second one P of many radio pulsars 
(see e.g. [151-453] and references therein). Two different classes of pulsars have been identified thanks 
to the observation of their rotational periods: the normal pulsars with rotational periods of the order 
of ~ 1s, and the so-called millisecond pulsars with rotational periods three orders of magnitude smaller. 
The first millisecond pulsar was discovered in 1982 with the help of the Arecibo radio telescope and 
since then more than 200 pulsars of this class have been observed. Nowadays, the fastest pulsar known 
until now, with a rotational period of 1.39595482 ms, is the object named PSR J1748-2446ad which 
was discovered in 2005 in the globular cluster Terzan 5 in the Sagittarius constellation. 
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Figure 12: Mass-radius relations obtained with different microscopic (solid lines) and phenomenological 
(broken lines) EoSs. The mass of the most heavy pulsar PSR J07404-6620 [330] observed until now is 
also shown, together with the constraints from the GW170817 event [8] and the mass-radius constraints 
of the NICER mission [432, 433]. 


6.2.5 Surface temperatures 


The effective surface temperatures of NSs can be determined from the detection of thermal photons 
from the stellar surface in isolated NSs (see Sec. 7.4) by fitting the observed spectra to blackbody ones. 
One should keep in mind, however, that NSs are not black bodies, because the hydrogen and helium (or 
even carbon) in their atmospheres modify the blackbody spectrum. Furthermore, the surface emission 
can be modified by the presence of strong magnetic fields. In fact, when realistic atmosphere models 
are used in the fit of the measured spectrum, surface temperatures are reduced. We note that any 
uncertainty in the determination of the surface temperature changes the corresponding luminosity Z 
of the star by a large factor according to the Stefan—Boltzmann law, L = 47 R?oT*. Consequently, it 
is not appropriate to use the surface temperature when comparing with observational data, but the 
luminosity instead. 


6.2.6 Gravitational redshift 


A source of very valuable information on NS structure is provided by the measurement of its gravitational 


redshift 
2GMN 12 
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which allows to constrain the M/R ratio and therefore the nuclear EoS. The interpretation of measured 
y-ray bursts, as gravitationally redshifted 511 keV e* annihilation lines from the surface of NSs, supports 
a NS redshift range of 0.2 < z < 0.5 with the highest concentration in the narrower range 0.25 < z < 
0.35 [454]. 


6.2.7 Quasiperiodic oscillations 


The observation and analysis of quasiperiodic x-ray oscillations (QPOs) [455] in x-ray binaries can 
provide very useful information to understand better the innermost regions of accretion disks as well as 
to put stringent constraints on the masses, radii and rotational periods of NSs. They may also serve as 
a unique proof of strong-field general relativity. QPOs measure the difference between the rotational 
frequency of the NS and the Keplerian frequency of the innermost stable orbit of matter elements in 
the accretion disk formed by the diffused material of the companion in orbital motion around the star. 
However, their theoretical interpretation is not simple and remains still controversial. 


6.2.8 Magnetic fields 


Since the pioneering work of Gold in 1968 [156] pulsars are generally believed to be rapidly rotating NSs 
with strong surface magnetic fields. The strength of the surface field could be of the order of 105 — 10? G 
in the case of millisecond pulsars, about 10!? G in normal pulsars, or even 10!4 — 10!° G in the so-called 
magnetars. The observation of the pulsar rotational period P and its first derivative P can be used to 
estimate the strength of the surface magnetic field of a pulsar, e.g. within the so-called magnetic dipole 
model [457] through the relation 

27? B? R6 sin? a 
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where B is the strength of the magnetic field at the pole, R is the radius of the NS, a is the angle 
between the magnetic and the rotation axis, and J is the moment of inertia of the star. 


6.2.9 Glitches 


Pulsars are observed to spin down gradually due to the transfer of their rotational energy to the emitted 
electromagnetic radiation. However, sudden jumps AQ of the rotational frequency Q have been observed 
in several pulsars followed by a slow partial relaxation that can last days, months or years. These jumps, 
mainly observed in relatively young radio pulsars, are known as glitches. The relative increase of the 
rotational frequency AQ/Q varies from ~ 107% to ~ 5 x 1076. The first glitches were detected in the 
Crab and Vela pulsars [158-460]. Nowadays we know more than 520 glitches in more than 180 pulsars. 


6.2.10 Timing noise 


The rotational stability is one of the most remarkable features of radio pulsars. Some of them, however, 
show slow irregular or quasi-regular variations of their pulses over time scales of months, years and 
longer times which have been called pulsar timing noise. These timing imperfections appear as random 
walks in the pulsar rotation (with relative variations of the rotational period < 10 !? — 1078), the 
spin-down rate or the pulse phase, which could stem from changes in the internal structure and/or the 
magnetosphere of NSs, although their nature is still uncertain and many hypotheses have been made, 
see e.g. Ref. [461]. 
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7 Equation of state and neutron star cooling 


We discuss in this chapter another physical phenomenon of NSs that is strongly related to the nuclear 
interactions and the EoS, namely the cooling properties of (isolated) NSs. Nuclear forces determine via 
the EoS not only the NS structure, but also the composition of NS cores, as well as reaction matrix 
elements, effective masses, and superfluid gaps of baryons, which are crucial for the computation of 
heat capacity and neutrino emission rates in NS cores. These two quantities govern the cooling of 
middle-aged (X 10° yr) isolated NSs, which can be confronted with astrophysical data on luminosity 
and surface temperatures for NSs of known or estimated age. 

The problem of NS cooling is a very rich field of physics that has been covered in several specialized 
reviews [162-472] and it has become impossible to give a complete overview. We try to single out 
here only the essential aspects related to the nuclear interactions and EoS. In particular, we focus 
on the cooling of isolated NSs and disregard the even more challenging problem of the cooling of 
accreting NSs in X-ray transients in quiescence (qXRT) [464, 473]. NSs may also be reheated by various 
other mechanisms, e.g. [474]. Furthermore, only purely nucleonic NSs are considered, disregarding 
scenarios comprising hyperons and quark matter, and also the important effects of magnetic fields are 
not considered in this brief chapter. 


7.1 Thermal evolution of a isolated neutron star 


We start by sketching the thermal evolution of a isolated NS (see [463, 464, 467] for a detailed review). A 
proto-NS is formed in a supernova event with a high temperature T ~ O(10" K). The proto-NS becomes 
a NS when it gets transparent to the neutrinos that are formed in its interior. During ~ 100 years, 
the crust stays hot due to its low thermal conductivity, while the core cools by emission of neutrinos. 
Therefore, the core and the crust cool independently and the evolution of the surface temperature 
reflects the thermal state of the crust and is sensitive to its physical properties (see, e.g. [175-478] and 
references therein). Then the thermal evolutions of core and crust couple. The cooling wave from the 
core reaches the surface and the whole NS cools by the emission of neutrinos, mainly from the core. 
This is the neutrino cooling stage. The evolution of the surface temperature depends mainly on the 
physical properties of the core. As the temperature in the interior of the NS continues to decrease, 
the neutrino luminosity becomes comparable to the photon luminosity. The NS then enters the photon 
cooling stage and the evolution of the internal temperature is governed by the emission of photons from 
the surface and is sensitive to the properties of the outer parts of the star. 

Currently, all isolated NSs for which the surface temperature has been measured thanks to X-ray 
observations are at least ~ 300 years old (Recently the ALAM radio telescope discovered a red blob in 
the remnant of SN1987 which has a high possibility to be a NS. If correct, this NS would be the youngest 
star observed with 33 yr [479, 480].) and they are all in the neutrino- and photon-cooling stages (see, 
e.g. Table 1 in [481]). Therefore they have an isothermal interior, i.e. the redshifted internal temperature 
T;(t) = T(t, rje* (1) is constant throughout the interior for densities larger than p ~ 10!? g/cm?. Here 
T (r,t) is the local temperature and ó(r) the metric function. 

The cooling evolution of a spherically symmetric NS can be characterized by two equations: the 
equation of thermal balance and the equation of thermal transport. For a spherical shell in the star, the 
change of thermal energy is caused by the neutrino emission, the heat flux passing through the surface, 
and the possible heating sources, for example, by converting magnetic or rotational energy into thermal 
energy [465]. This is expressed by the thermal balance equation 
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where Q, stands for all relevant sources, Q, is the neutrino emissivity, C, the specific heat 

capacity, and m is the gravitational mass enclosed inside the sphere of radius r. The local luminosity 
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L is defined as the non-neutrino heat flux passing through the surface of the spherical shell. Since the 
heat flux is transported through thermal conduction, one can write the equation of thermal transport 


as 
L 2Gm o 
— = — L—— —ó fo) 
7772 KA/1 Ee 5, € T), (65) 


where & is the thermal conductivity. Eqs. (64) and (65) are partial differential equations for luminosity 
L(r,t) and temperature T (r,t), and all quantities Qa, Qu, Cy, x depend as well on the radial coordinate r 
and time t. On the stellar surface, luminosity and temperature can be connected by 


L (t) = L(R,t) = gend RT), (66) 


where L. and T, are photon luminosity and effective surface temperature, respectively, and dsg is the 
Stefan-Boltzmann constant. 

Once the interior of the star has become isothermal (T;) after about 10? years, isolated by a heat- 
blanketing atmosphere with a effective surface temperature 7;, the above equations simplify to the heat 
equation 

dT; 
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where the relation 7,(T;) is determined by the model of the heat-blanketing envelope, see Sec. 7.3. C(T;) 
is the total specific heat integrated over the whole star and the quantities L>°, D are respectively the 
redshifted total neutrino luminosity and the surface photon luminosity detected by a distant observer. 


IAN), (67) 


Observable quantities for astrophysical observation at large distance are apparent radius Rx, effec- 


tive surface temperature 7°, and apparent photon and neutrino luminosities L7, L> [482], 


Ro = RA —2GM/R, (68) 
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with M = m(R) the gravitational mass of the star of radius R. Hence, by solving Eqs. (64,65) or 
Eq. (67), one can obtain the dependence of the apparent photon luminosity L (or of T?°) on the 
stellar age t. This is the main goal of the cooling theory. 


7.2 Models of the crust 


The crust of an isolated NS born in a core-collapse supernova explosion is formed during the cooling of 
the very hot and fully fluid proto-NS [483, 484]. In the initial state, the composition of the outer layer 
corresponds to nuclear equilibrium, because T > 10!°K and therefore nuclear reactions are sufficiently 
rapid. The hot plasma crystallizes while cooling, and a standard assumption is that during the process 
of cooling and crystallization the plasma keeps the nuclear equilibrium. Consequently, when matter 
becomes strongly degenerate, the structure and EoS of the crust can be well approximated by cold 
catalyzed matter, the ground state of the matter at T = 0. The structure of this crust may be 
calculated to various degrees of sophistication, as reviewed in more detail in Sec. 4. In practice, the 
contribution of the crust to the neutrino cooling is usually negligible. 
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Figure 13: Luminosity (left panel) and temperature (right panel) vs age of 55 isolated NSs. 
Adapted from Ref. [481]. Numbers refer to Table 1 of that reference. See also the online database 
http://www.ioffe.ru/astro/NSG/thermal/cooldat.html (where object 1 has been withdrawn). 


7.3 Heat-blanketing envelope 


Within a few hundred years the redshifted temperature inside a newly born NS for densities larger than 
py ~ 1019 g/ cm? becomes uniform due to the high thermal conductivity. However, in the NS atmosphere 
the heat transport is dominated by the photons and in between there exists a thin layer, which has a 
low thermal conductivity, since the electrons are not highly degenerate and the large density strongly 
prevents photon transport. This results in high temperature gradients in the envelope, that is a few 
hundred meters thick. 

Therefore, a variety of models are devoted solely to the precise modeling of the heat-blanketing 
envelope [485-488], in the plane-parallel and stationary approximation, resulting in a relation between 
the surface temperature T, and the internal temperature 7; at oy. Some models consider various 
compositions for the envelope and, in particular, different abundances of light elements such as hydrogen, 
helium, carbon, resulting from the accretion of matter [486] or binary mixtures of these [479]. The higher 
the abundances of light elements, the higher the surface temperature for a given T;. As these abundances 
are unknown for a given NS, in the simulations one usually considers two limiting cases corresponding to 
the absence of light elements (non-accreted envelope) and a maximum amount of them (fully-accreted 
envelope), in order to estimate the possible range of predictions. 


7.4 Luminosity data 


The observable luminosity of isolated NSs is basically the photon luminosity, since the neutrino lu- 
minosity is far too weak to be detected directly. However, due to the unknown composition of the 
atmosphere just discussed, the unknown distance, or the interstellar absorption, it is difficult to obtain 
accurate results of photon luminosities. Even so, these data are still essential to impose constraints 
on the theoretical study of NS cooling. Of course, also the ages of neutron stars need to be known. 
They can be extracted from the spin period and its time derivative, which are the characteristic ages 
t = P/2P [489], or from kinetic properties of the stars (proper motion for instance), which are the 
kinetic ages. The kinetic ages are favoured where possible, and characteristic ages are treated as upper 
limit in most cases. Compilations of these data have been given in [464, 467-469, 473, 481, 490, 491] 
and are continuously updated. In Fig. 13 we reproduce the most recent set [481] of the bolometric 
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luminosities L5° and surface temperatures T?° of 55 NSs. It has too be stressed and can (not) be seen 


in the figure that weak sources are simply too faint to be observed and therefore these data are subject 
to an important and difficult to quantify selection effect, which strongly hampers the interpretation of 
theoretical cooling models that rely on unbiased data. 


7.5 Microphysical ingredients 
7.5.1 Nucleon effective masses 


The most prominent input information regarding single-particle propagation are the in-medium nucleon 
(n, p) effective masses m* = kp/[de(k)/dk], related to the density of states m*kr/7?. As most other 
microphysics input depends on them, they affect in several instances the computation of NS cooling. 
Note that these are already in-medium modified quantities, whereas other ingredients (matrix elements, 
cross sections, ...) are still often approximated by vacuum results, as will be discussed in the following. 

In principle, they should therefore be computed in a consistent manner together with the EoS. 
However, since other ingredients of the cooling calculations (matter composition, matrix elements, 
pairing properties) are far more uncertain and grossly affect the final results, this is not always done 
and bare nucleon masses are simply used, or different approximate but inconsistent values are used in 
the various components of the cooling simulations. In any case, in several theoretical approaches for the 
EoS the effective masses can be calculated straightforwardly [60, 492, 493], even at finite temperature 
[494-496]. We refer to the specialized review [497]. 


7.5.2 Specific heat and thermal conductivity 


The main contribution to the overall specific heat C(T;) in Eq. (67) comes from the NS core, i.e. from 
the neutrons, protons, and electrons. If the core is non-superfluid, most of the specific heat comes from 
the neutrons, as indicated by the simplest result for a Fermi gas at density n and ‘low’ temperature 
T «ep, 


p (72) 


For the less important crust contributions there are several refinements [169]. If the neutrons and protons 
are superfluid, then below their superfluid critical temperature 7; their values become exponentially 
suppressed [162, 498] and C is dominated by the electrons [488]. 

The local thermal conductivity & in Eq. (65) is a sum of contributions of the different constituents 
of matter. Here we only illustrate qualitatively the dominant neutron contribution in the core [469, 
471, 488, 499—501], neglecting the coupling to protons, 


1 57? m?n 


~ Lon 73 
KE F192 m*4o ' (73) 


where o is the effective cross section at the root of the phenomenon. In-medium effects [60, 501-505] 
provide an effective mass m* < m, but also modify the effective cross section in the opposite way, 
o" > o. The overall effect is an enhancement of the conductivity [501, 503, 506-508], which remains 
dominant compared to the lepton conductivity [488, 509, 510] for temperatures T < 10% K. However, 
the theoretical results of different groups differ by even an order of magnitude, see the discussions in 
[500, 501]. As a further complication, superfluidity affects the thermal conductivity in a highly nontrivial 
way and can either enhance or reduce it [471, 500]. 

The quantitative understanding of the thermal conductivity is only important for young NSs, t S 
O(10? yr), where thermal relaxation is not yet reached in the interior. For example, a reduction of 
the thermal conductivity was assumed in the medium-modified one-pion-exchange (MOPE) model of 
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Figure 14: Elementary diagrams for neutrino emission. 


[511—514] in order to explain the apparent accelerated cooling of the Cas A NS [515-520] with an age 
of about 330 years. 


7.5.3 Neutrino emissivity 


In this section we briefly recall the main neutrino emission mechanisms in the NS core, which dominate 
the total neutrino luminosity L>°, and the relevance of the nucleon effective masses and other in-medium 
effects, following closely the detailed treatment given in Ref. [463] and updates in [469, 471]. Only the 
rates for the non-superfluid scenarios will be given in this section, for which the dependence on the 
effective masses is via the general factor [60] 


1/3 «Ni *\ J 
me (2) (aa) Ge) m 
Po Mn Mp 


In the presence of superfluidity the dependence becomes highly nontrivial and requires detailed calcu- 


lations [463]. In the following all emissivities Q, are given in units of ergem”*%s”?, 


Direct Urca 

In the absence of pairing three main mechanisms are usually taken into account, which are the direct 
Urca (DU), the modified Urca (MU), and the NN bremsstrahlung (BS) processes, as sketched in Fig. 14. 
By far the most efficient mechanism of NS cooling is the DU process, Fig. 14(a), which is in fact the 
in-medium neutron beta decay followed by its inverse reaction [216, 521]: 


n-— pde uw and pte >n+. (75) 
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However, the energy and momentum conservation imposes a density threshold on this process. The 
result for npe NS matter is given by 


QU) a 4.0 x 10M TOP + KO — KØ), (76) 


where Ty is the temperature in units of 10? K. If muons are present, then the corresponding DU process 
may also become possible, in which case the neutrino emissivity is increased by a factor of 2. Thus 
the DU process appears in a step-like manner, once the proton fraction exceeds a value of 13-14 % 
(slightly dependent on the muon fraction at onset). As the proton fraction x is determined by the 
nuclear symmetry energy, 


(37? px) "3 = pe = us — Mp SAL — 20) Esym(p) ; (77) 


EoSs with a large enough symmetry energy feature DU cooling. This is the case for many microscopic 
EoSs, see Fig. 2 for illustration. 


Modified Urca 

If the DU process is kinematically forbidden or strongly reduced, various less efficient neutrino 
processes may be operating in the NS core. While the former is a pure weak reaction, those processes 
are driven by strong interactions, and are consequently affected by much greater theoretical uncertainties 
that will be discussed in the following. The two main ones with an emissivity Q, « T® are the MU 
processes, Fig. 14(b-e), 


N+n>N+p+e tv, and N+p+e 3 N+n+v, (78) 


where N = n,p is a spectator nucleon that ensures momentum conservation. Since five degenerate 
fermions are involved instead of three, the efficiency is significantly reduced compared to the DU process. 
The emissivities of the MU processes in the neutron and proton branches, employing for the long-range 
part of the in-medium NN interaction (gray block in Fig. 14) the dominant one-pion-exchange (OPE) 
contribution and for the short-range part an effective contribution in the framework of Landau theory, 
are given by [463, 469, 522, 523] 


QUI) = 8.1 x 10 Mai TB, ; (79) 
QM?) 2: 8.1 x 10% MT 098p (1 — KEARE) OKP +10 — Ki), (80) 


where the factors a, a, take into account the momentum-transfer dependence of the squared reaction 
matrix element under the Born approximation, and Bn, 8p include the non-Born corrections due to 
NN interaction effects, which are not described by the OPE [463]. The currently adopted values are 
Ap = An = 1.13 and p = n = 0.68. The main difference between the proton branch and the neutron 
branch is the threshold character, although usually irrelevant. If muons are present in the dense NS 
matter, the equivalent MU processes become also possible, and accordingly several modifications should 
be included in Eqs. (79,80), as discussed in Ref. [463]. 

As stated above, these MU rates have been derived using a very simple approximation (free OPE + 
Landau parameters) for the in-medium NN interaction in the relevant diagrams Fig. 14(b,c), and are 
consequently density independent apart from the effective-mass corrections and kinematic factors. This 
has been attempted to improve in several works. 

In particular, in the medium-modified OPE (MOPE) model of [524-527] one asserts that the emis- 
sivity receives a strong and dominant density-dependent correction due to the softening of the pion 
mode (reduction of the in-medium pion mass). Then at p > po, diagrams of type Fig. 14(e) dominate 
(not included in the OPE result), which describe in-medium conversion of a virtual charged pion to a 
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neutral pion with the emission of a real lepton pair. The modification factor for the medium-modified 
Urca reaction with respect to the free OPE result is (511-513, 528-530] 


pMMETS " (£y ea) (81) 


where T(p) © 1/[1 + 1.6(p/po)'P] is an estimate of the effect when replacing the bare TNN vertex by 
the dressed one, computed from Landau-Migdal parameters, and 

&^ = min [—D*(w = 0, EI (82) 
is the effective pion gap in the medium, computed from the dressed pion propagator D,(w,k). The 
density dependence of w is model dependent and also related to a possible pion condensation at high 
density. With typical parameters, one obtains enhancement by a factor of ~ 3 at po and up to 5000 at 
3po. Note, however, that this enhancement strongly depends on the uncertain values of the pion gap 
and the vertex correction at large density. 

On the contrary, when replacing the OPE approximation by the in-medium T-matrix [531], a re- 
duction to about 1/4 of the OPE result was found. 

Recently [532], a further ‘kinematic’ in-medium effect was explored, namely that for densities above 
the DU threshold the propagator of the virtual nucleon after emission of the lepton pair in Fig. 14(b) 
might develop a pole that is unaccounted for in the standard treatment leading to Eqs. (79,80). This 
correction of the MU rates leads to a significant enhancement by a factor of several, in particular close 
to the DU threshold p = ppy, and thus provides a more smooth switch-on of the DU process. 

The problem of an accurate computation of the MU rates can currently be considered unresolved, 
as no formalism incorporating consistently (also together with the computation of the EoS) the various 
in-medium effects has yet been attempted. It seems, however, that the simple expressions Eqs. (79,80) 
constitute underestimates of the real MU rates. 


Bremsstrahlung 
In the absence of the DU process, the standard neutrino luminosity of npe matter is determined not 
only by the MU processes but also by the BS processes in NN collisions, also described by Fig. 14(b-e), 


N+NSN+N+v+7, (83) 


with N a nucleon and v, v an (anti)neutrino of any flavor. These reactions proceed via weak neutral 
currents and produce neutrino pairs of any flavor [522, 523]. Contrary to the MU, an elementary 
act of the NN BS does not change the composition of matter. The BS has accordingly no thresholds 
associated with momentum conservation and operates at any density in uniform matter. In analogy with 
the MU process, the emissivities depend on the employed model of NN interactions and in OPE+Landau 
approximation are [463] 


QU?) à, 2.3 x OMT Ona Ban ën d IT , (84) 
QEP 24 4.5 x 107 MT Deen Ban > 2 
QPPP) = 2.3 x 10% Mad TS o, pp - ES 


As for the MU rates, the dimensionless factors ayy come from the estimates of the squared matrix 
elements at p = po: Ann = 0.59, Qnp = 1.06, App = 0.11. The correction factors nn are taken as 
Ban = 0.56, Bnp = 0.66, Bp, = 0.70. All three processes are of comparable intensity, with QPPP) < 
Qr) < Or Bon). The simple estimates Eqs. (84,85,86) are about 50 times weaker than the OPE MU 
rates Eqs. (79,80). However, also in this case the in-medium effects might completely change the picture: 
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As for the MU rates, calculations show that the use of the realistic T-matrix instead of the OPE 
leads to the suppression of the neutrino emissivity approximately by a factor of four [533-537], although 
[531] found an even larger decrease by a factor of 10-20 in the medium. Several medium effects to 
second order in the low-momentum universal potential Vowk were evaluated in [538], yielding an overall 
reduction of the emissivity by a density-dependent factor of about 2-4, compared to the standard factor 
Canon 0.35. In the MOPE model the in-medium correction for BS is [511-513, 529, 530] 


MBS Su p 4/3 [T (p) /T(oo)]9 
el) me (87) 


Contrary to the T-matrix corrections, significant enhancement at p > po up to a factor of 8489 ~ 100 
is obtained, depending on the model adopted for the pion propagator. 

We summarize that the strong nuclear interaction influences directly several ingredients of the MU 
and BS neutrino emission rates, and very large in-medium effects have been claimed in some publica- 
tions. However, conflicting results of both enhancement and suppression are reported depending on the 
specific medium effect considered, and no generally accepted view has yet emerged. However, of even 
greater importance for the cooling theory are the superfluid properties of nuclear matter, as discussed 
in the following. 


7.6 Effects of superfluidity 


The neutrino cooling processes in NSs are dramatically affected by the neutron and proton superfluidity 
[462, 463, 472, 539-544], and the knowledge of the pairing gaps A in the 150 and 3PF2 channels of the 
NN interaction in beta-stable matter is essential for a correct description of the thermal evolution of 
a NS. These superfluids are produced by the pp and nn Cooper pairs formation due to the attractive 
part of the NN potential, and are characterized by a critical temperature Te = 0.567A for isotropic 
gaps. The pairing gaps reduce the phase space available for the cooling reactions and also modify the 
neutrino-emission matrix elements by the effect of anomalous propagators. The latter feature is nearly 
always disregarded, leading to universal superfluid modification factors which do not depend on details 
of the strong-interaction matrix elements. The occurrence of pairing leads to two relevant effects in NS 
cooling, namely 

(i) A strong reduction by several orders of magnitude when T < T; of the emissivity of the neutrino 
processes the paired component is involved in, with a corresponding reduction of the specific heat of 
that component (Sec. 7.5.2). For example, proton superfluidity in the core of a NS suppresses both 
Urca processes (DU and MU) but does not affect the neutron-neutron BS. The general expressions 
cannot be obtained analytically, but fit formulas are available (163, 545]. For illustration we list the 
suppression factors of the emissivities R = Qa/Qo in the case of p1S0 pairing in the low-temperature 
limit y = A/T > 1 for the various neutrino reactions [463, 523, 545, 546]: 


Rpy = 7.19 x 104447, (88) 
Rmn © 0.166 x 10 "wie", Ru, ™ 144 x 10° ve”, Ge 
Hgap Si 0516ve*”, Rpgpp 52 0.216 2e 7" , (90) 


At intermediate temperatures, the suppressions follow power laws, though. Note that in the absence of 
neutron pairing, the Bnn reaction becomes the dominant unsuppressed cooling mechanism. 
(ii) Onset of the *Cooper pair breaking and formation" (PBF) process with associated neutrino- 
antineutrino pair emission, 
N3N+v+9, (91) 


which is in fact a neutral-weak-current variant of the DU process Fig. 14(a), and kinematically forbidden 
in normal matter. However in superfluid matter it becomes unblocked [463, 547, 548] and the emissivity 
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can be written as eg 
m 
QUPD e 3.5 x 1091 NE Tay Fw (y) (92) 
MN 
with y = A(T)/T. For both 150 and 3P2 pairings, the emission of neutrinos by the PBF process 
can occur through axial and vector channels. The computation of the factor ay = ak + aj, from 
the elementary superfluid weak interaction matrix element and its renormalization due to the strong 
interaction is very intricate and has been tackled during the last decade [463, 549-564], with the current 
state of results for p1S0 and n3P2 pairing: 


ap X 0,0: + SE e hn” 94/2 x 0.8, (93) 


where gy = 1 and ga = 1.26 are the nucleon weak vector and axial-vector coupling constants. Thus, 
the PBF process for singlet pairing is a relativistic effect (vp = kp/m* < 1) and much weaker than 
the other cooling processes. However, the PBF contribution for triplet pairing is important and even 
dominant for certain temperature ranges, as determined by the universal control function 


2 
oo r? + y? 
Eye a |). (94) 
0 1 even? 


Qualitatively this process starts once the temperature decreases to T, of a given type of baryons, becomes 
maximally efficient when T z 0.8 T, ~ 0.45A, and then is exponentially suppressed for T < T, [463]. 
Thus during the cooling evolution and in regions of the star where locally T(r,t) = 0.4A(r), the n3P2 
PBF process might be the dominant cooling reaction. 


7.7 Nuclear pairing gaps 


As evidenced in the previous section, the most important nuclear physics input for the superfluid 
properties of stellar matter are the pairing gaps A;(p) in the different partial waves. Usually the most 
important ones are the p1S0 and n3P2 pairing channels, while the n1S0 (crust only) gap is much less 
relevant for the cooling [565], and the p3P2 gap is normally disregarded due to its uncertain properties 
at extreme densities. 

Microscopically the gaps are determined by the irreducible NN interaction in the different partial 
waves, apart from the s.p. properties (EoS). In the simplest BCS approximation, and detailing the more 
general case of pairing in the coupled 3PF2 channel, the pairing gaps are computed by solving the 
(angle-averaged) gap equation [292, 566-570] for the two-component L = 1,3 gap function, 


Fe (k) = - 1 "akk am (rå (k, k’) = (k^) (95) 


with 
E(k}? = [e(k) — nu + Ar (1)? + As(k)?, (96) 
while fixing the (neutron or proton) density, 
_ ke oyo ll, bin 
IS ER) | On 


Here e(k) = k?/2m + U(k) is the s.p. energy, u = e(kr) is the chemical potential determined self- 
consistently from Eqs. (95-97), and 


Vp (kk) = / dr? jp (Kr) VIS (0) jo (kr) (98) 
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Figure 15: Theoretical predictions for the 150 gap in pure neutron matter [502, 572-575, 577, 579-582, 
584-587, 590]. Curves represent traditional approaches and markers denote Monte Carlo calculations 
of different kind. 


are the relevant matrix elements (J = 1 and S = 1; L, L’ = 1,3 for the 3PF2 channel, S = 0; L, I/ =0 
for the 150 channel) of the bare NN potential V. The relation between (angle-averaged) pairing gap 
at zero temperature A = y A%(kp) + A3(kp) obtained in this way and the critical temperature of 
superfluidity is then T, = 0.5674. (If no angle average is performed, the prefactor varies slightly, see, 
e.g. (162, 463], but the angle-average procedure is usually an excellent approximation [567, 571]). 


However, in-medium effects might strongly modify these BCS results, as both the s.p. energy e(k) 
(effective mass m* and quasiparticle strength Z) and the interaction kernel V itself are affected by TBF 
(Sec. 2.1.2) and polarization corrections. It turns out that in the p1S0 channel all these corrections lead 
to a suppression of both magnitude and density domain of the BCS gap [81, 502, 572-590], see Fig. 15 
for a compilation of previous results and, e.g. [172, 539, 586, 591] for further reviews. The results are 
widespread, which certifies the extreme sensitivity of the gap to details of interaction, approximations, 
and theoretical approach. 


The situation is even worse for the gap in the n3P2 channel, which already on the BCS level 
depends on the NN potential (292, 568, 570, 588, 592, 593], as at high density there is no constraint 
by the NN phase shifts. Furthermore TBF act generally attractive in this channel, but effective mass 
and quasiparticle strength reduce the gap and polarization effects on V might be of both signs, in 
particular in asymmetric beta-stable matter [571, 578, 585, 588, 594-597]. Note that most theoretical 
investigations consider only pure neutron matter. Thus, due to the high-density nature of this pairing, 
the various medium effects might be very strong and competing, and there is still no reliable quantitative 
theoretical prediction for this gap. 


7.8 Cooling simulation scenarios 


As detailed in the preceding sections, the sophistication of cooling simulations has evolved with time 
and we briefly report here the main developments. 
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Figure 16: Example cooling diagrams for various cooling scenarios: (a) Minimal cooling scenario 
with the APR EoS [95] (Mpu = 2.03 M5) and heavy-element envelopes. For the cooling curves of 
M = 1.4 Mọ, three choices of the n3P2 gaps are considered (‘0’, ‘1’, ‘2’; the size of the n3P2 gap 
increases with the number). Adapted from Ref. [490]. (b) Including different DU processes with the 
Gan EoS [598] (Mpu = 1.38 Mo) inside a 1.4 Mo NS, except the kaon case whose mass is 1.8 Mo. 
Adapted from Ref. [599]. (c) Nuclear medium scenario with the HHJ (parameterized softened APR) 
EoS [600] (Mpy = 1.84 Mo) for several choices of the NS mass. The thin red curves indicate the cooling 
without pairing, while the thick black curves employ both neutron and proton gaps, but the n3P2 gap 
suppressed by a factor 10. Adapted from Ref. [511]. 


7.8.1 Minimal cooling 


The occurrence of superfluidity in NS matter is generally accepted and its effects have to be considered in 
the NS cooling simulations. A proposed extension of the standard cooling scenario without superfluidity, 
the so-called ‘minimal cooling scenario’ [490, 555], includes the neutrino emission from the PBF processes 
as the fastest cooling reaction (see Sec. 7.6) and assumes that other enhanced neutrino emissions (various 
DU processes) are not allowed. This assumption is justified if stellar matter has a too small proton 
fraction (X; 13%), and furthermore no exotic matter is present in the core. Within minimal cooling 
scenarios, the MU processes are suppressed by neutron and proton pairing gaps, and the n3P2 PBF 
process dominates cooling in the entire neutrino cooling era. This remains valid even with different 
models of pairing gaps [190]. The p1S0 PBF reaction is also present in the core of NSs, however, its 
effect is practically negligible [490, 555, 601]. 


A detailed discussion of the minimal cooling is given in Refs. [490, 555]. An example of cooling 
curves with PBF cooling [555] is shown in panel (a) of Fig. 16. With PBF cooling NSs have a lower 
luminosity compared to the one without superfluidity, and this depends on the sizes of the gaps. In 
addition, the cooling curves have a weak dependence on the stellar mass. This results in a difficulty 
to describe some young cold NSs and old hot NSs. Accordingly, for the minimal cooling model to be 
consistent with data, it has to be fine-tuned to assure the most efficient cooling [465, 514, 555, 602], 
although fine-tuning the n3P2 gap [603] allows also to reproduce the speculated accelerated cooling of 
the supernova remnant Cas A [515-520]. 
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Figure 17: Cooling diagram for the BHF V18 EoS [191] with p1S0 BCS gap and no n3P2 pairing, 
for different NS masses M/M¿ = 1.0,1.1,...,2.3 (decreasing black curves). The dashed green curve 
(hardly visible) marks the NS mass Mpu + 0.01 Mọ = 1.02 Mọ at which the DU process has just set 
in, the dash-dotted green curve marks the NS mass Mjg9 = 1.92 Mo for which the p1S0 gap vanishes 
in the center of the star, and the dotted green curve corresponds to Mmax = 2.34 Mo. The black 
curves are obtained with a Fe atmosphere and the shaded areas cover the same results obtained with a 
light-elements (7 = 107) atmosphere. The data points are from [481]. Adapted from [601]. 


7.8.2 Importance of DU cooling 


The efficient DU process seems to be required in a consistent cooling model in order to reproduce 
both relatively hot (XMMU J1731-347, no. 10 in Fig. 13) and relatively cold (PSR J0205+6449, 
PSR B2334+61, PSR J0007+7303, PSR B1727-47; nos. 14,33,46,48) cooling objects in the current 
database. Apart from the nucleonic DU process, possible exotic matter in the NS core, e.g. hyperons 
[233, 604, 605], pion or kaon condensation [216, 511], and deconfined quarks [529, 606, 607], could also 
provide DU processes and result in a fast cooling, as illustrated in panel (b) of Fig. 16. 

The possibility of a high enough proton fraction to allow the nucleonic DU process was realized early 
in schematic models [216, 521, 608]. In fact many modern microscopic nuclear EoSs do reach easily the 
required proton fractions for the DU process [42, 95, 121, 190, 599, 609] and are able to describe current 
cooling observations of isolated NSs [565, 601, 610], as well as the reheated cooling of the accreting 
NSs in X-ray transients in quiescence [478, 611, 612], provided the DU process is dampened by the 
presence of nuclear pairing throughout a sufficiently large (but not complete) set of the NS population. 
In this framework a successful consistent modeling of all cooling data (in particular the unusually hot 
XMMU J1731-347 [613]) requires a p1S0 gap extended over a large enough density/mass range and a 
very small or vanishing n3P2 gap [473, 481, 511-513, 565, 601, 610, 611, 614, 615]. This is opposite to 
the minimal cooling scenario, where n3P2 PBF cooling is required as the only fast cooling reaction. 
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Fig. 17 shows an example obtained with the BHF V18 EoS (Sec. 3.3) [601], featuring a threshold 
mass of DU cooling Mpy = 1.01 Mo, so that practically all cooling curves involve nucleonic DU cooling. 
The main effect of superfluidity on NSs with M > Mopv is the strong quenching of the DU process, and 
thus a substantial reduction of the total neutrino emissivity. Hence those stars have a higher luminosity 
than in the non-superfluid case. On the other hand, if M > Mıso = 1.92 Mo (corresponding to a NS 
with the p1S0 gap vanishing in the center), the complete blocking of the DU process disappears and 
the star cools very rapidly again. One observes results in line with these features in the figure, namely 
between Mpy and Mjgo there is a smooth dependence of the luminosity on the NS mass for a given 
age. All objects in the database can be explained by choosing a proper atmosphere model in each case. 


7.8.3 Medium-modified Urca processes 


As discussed in Sec. 7.5.3, the standard cooling rates disregard various in-medium effects. The authors 
of the “nuclear medium cooling” approach [511, 526-528, 532] advocate very strong modifications of 
the different cooling processes, in particular a strong density-dependent enhancement of the MU and 
BS reactions, Eqs. (81,87). This causes a smooth crossover from the standard to the enhanced cooling 
scenario for increasing NS masses and can describe reasonably well the young cold NSs without the 
inclusion of DU processes and relevant gaps (511, 513, 514]. Panel (c) of Fig. 16 shows that cooling 
curves of 1.0 Mo (solid curve) and 1.4 Mo (dashed curve) NSs are widely separated compared to minimal 
cooling where the curves overlap. But, superfluidity is required for describing hot NSs. The slow cooling 
and intermediate cooling objects (for instance object 2, RX J0822-43) cannot be reproduced without 
gaps (thin curves). However, the inclusion of the n3P2 gap results in a too efficient PBF process and 
a too fast cooling, which prevents to explain at least several of the cooling data. One requires thus a 
strong suppression of the n3P2 gap [511, 513, 514, 604, 614], just as in the models involving DU cooling. 


7.9 Epilogue 


Accurate quantitative modeling of NS cooling with account of all possible effects, even on the purely 
nucleonic level, is very complicated and has so far only been achieved in rather simplified ways. We 
reiterate here only two of the most critical issues: 

(1) The elementary cooling rates, heat capacities, thermal conductivities, etc., have been computed 
based on very simplified in-medium NN interactions (e.g. OPE with severe kinematic averages). Going 
beyond this approximation in one way or another (T or G matrix, in-medium pion propagation, kine- 
matic corrections, relativistic corrections, ...) seems to indicate very important in-medium corrections 
(potentially of several orders of magnitude), but no consistent approach (using an in-medium interac- 
tion computed consistently based on the same NN interaction as for the nuclear EoS) has yet been 
presented to our knowledge. Therefore most current theoretical results may in general be considered 
only order-of-magnitude estimates. 

(ii) An even more important regulator of cooling rates are the neutron and proton pairing gaps in 
NS matter. While there is some general agreement on a substantial suppression of the p1S0 gap, the 
high-density n3P2 gap is currently theoretically not under control. 

The main issue in cooling simulations is whether or not the proton fraction in the stellar matter is 
high enough to allow DU cooling. The former is the case for many modern microscopic EoSs. Then 
cooling is determined practically only by properties of the DU process. Evidence has been put forward 
for a theoretically-supported nucleonic scenario in which the DU process is active for high-mass stars, 
and damped by p150 pairing for intermediate-mass stars, while n3P2 pairing is not allowed. However, 
alternative scenarios are still possible and can not generally be excluded by the current limited amount 
of cooling data. 

Going beyond the nucleonic level, the cooling can be affected by the presence of hyperons (see Sec. 5 
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and, e.g. [233, 468, 604] and references therein) or deconfined quarks ([529, 565, 606, 607] and references 
therein), by pion or kaon condensation (see [216, 463, 511] for review and references), by emission of 
axions (e.g. [616] and references therein), and other effects. These scenarios are currently even more 
uncertain and unexplored than the purely nucleonic one. 


8 Conclusions 


In this review, we have presented the current status of the EoS for compact objects, and discussed 
the different theoretical approaches, both microscopic and phenomenological, for constructing the EoS. 
Only nucleonic and hyperonic degrees of freedom were taken into account. In spite of the tremendous 
progress in the last few years regarding the actual constraints on the EoS, coming from both nuclear 
physics experiments and astrophysical observations, especially thanks to the recent GW detection, it is 
not possible yet to seriously exclude one or more scenarios depicted by the different theoretical models. 
Recently a lot of attention has been devoted to the results on NS mass and radius from the NICER 
mission, as well as the tidal polarizability measured for the GW170817 event, along with the observation 
of the large mass of the millisecond pulsar MSP J0740+6620, M = 2.141938 Mo. These constraints are 
at the moment the more stringent for selecting the EoS, whereas those coming from laboratory nuclear 
experiments are less strict since often model dependent. 

We expect that in the near future a new generation of telescopes and projects such as eXTP [617] 
and advanced GW detectors such as LIGO and Virgo can provide more and more precise data that 
can significantly contribute to probe the internal structure of compact objects, thus disentangling the 
different theoretical models. Besides that, we have shown that the cooling theory of NSs is a rich field of 
physics with a great potential to constrain the high-density stellar EoS. Here the main issue is whether 
the EoS is stiff enough to allow a high proton fraction and the rapid DU cooling process that dominates 
all other mechanisms. Currently theoretical models are strongly hampered by uncertainties mainly 
related to a correct treatment of the in-medium interactions, so that theoretical predictions must be 
considered much less reliable than those related to the ‘ordinary’ EoS at subnuclear densities. 

Despite all efforts, we remark that there are still several open issues in the EoS modelling, which 
need to be addressed in the near future. These include: 


e The role of TBF in the nuclear-matter EoS if only nucleonic degrees of freedom are considered; 


e The solution of the “hyperon puzzle” in microscopic approaches and the more general question 
about the treatment of non-nucleonic degrees of freedom; 


e The conditions for a phase transition to quark matter to take place; 
e The importance of thermal effects for the construction of the EoS; 


e The cooling mechanisms, in particular pairing, consistent with the EoS. 


Answers to those questions will surely improve our understanding of the various fundamental phe- 
nomena related to the astrophysics of NSs, their internal constitution, the explosion mechanisms of 
core-collapse supernovae, the mass threshold for BH formation, the dynamics of NS binary mergers, 
and the nucleosynthesis of heavy elements. Therefore the current theoretical uncertainties will require 
significant efforts to be undertaken in these directions in the near future. 
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